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Abstract—Antenna current optimization is a tool that offers 
many possibilities in antenna technology. Optimal currents are 
determined in the antenna design region and used for physi- 
cal understanding, as a priori estimates of the possibilities to 
design antennas, physical bounds and as figures of merits for 
antenna designs. Antenna current optimization is particularly 
useful for small antennas and antennas that are constrained 
by their electrical size. The initial non-convex antenna design 
optimization problem is reformulated as a convex optimization 
problem expressed in the currents on the antenna. This convex 
optimization problem is solved efficiently at a computational cost 
comparable to a Method of Moments (MoM) solution of the same 
geometry. 

In this paper a tutorial description of antenna current opti- 
mization is presented. Stored energies and their relation to the 
impedance matrix in MoM is reviewed. The convex optimization 
problems are solved using MATLAB and CVX. MoM data is 
included together with MATLAB and CVX codes to optimize the 
antenna current for strip dipoles and planar rectangles. Codes 
and numerical results for maximization of the gain to Q-factor 
quotient and minimization of the Q-factor for prescribed radiated 
fields are provided. 


Index Terms—Stored energy, Antenna Q, Convex Optimiza- 
tion, Physical bounds, Q-factor, Antenna theory 


I. INTRODUCTION 


Antenna design can be considered as an art to shape and 
choose the material to produce a desired current distribution 
on the antenna structure. Antenna current optimization is a 
preliminary step where the current distribution is determined 
for optimal performance with respect to some parameters. 
This step splits the antenna synthesis process into two less 
complex tasks; the first one is to determine the optimal 
current and the second is to determine an antenna structure 
that provides similar performance to the optimal current. The 
current distribution serves as a guideline to antenna design 
but it is most useful for determining an upper bound on the 
antenna performance, i.e., a physical bound or fundamental 
limitation. 

Optimization is common in antenna design to augment 
existing structures and to construct new designs [62]. Meta- 
heuristic algorithms, such as genetic algorithms [44], particle 
swarm and gradient-based algorithms dominate the antenna 
optimization field due to the inherent complexity of antenna 


design problems. Optimization of the current density on the 
antenna is inherently different and can often be formulated 
as a convex optimization problem [31]. The formulation in 
convex form is advantageous because it covers a broad range 
of different problems by combining constraints. There are also 
many efficient solvers for convex optimization problems and 
these solvers can provide error estimates [6, 19]. Antenna 
optimization parameters can be combined as quadratic forms, 
such as stored energy and radiated power; linear forms, such 
as near- and far-fields and induced currents; and norms to 
formulate convex optimization problems relevant for a specific 
antenna problem [31]. 

One of the most challenging computational tasks in antenna 
current optimization is the evaluation of the stored energy [29, 
30, 70]. Fortunately, the matrices used to compute the stored 
energy are, in principle, already implemented in many Method 
of Moments (MoM) solvers. What is needed in the perfect 
electrical conductor (PEC) case is to separate the electric and 
magnetic parts of the impedance matrix from the Electric Field 
Integral Equation (EFIE) and to add a non-singular part. This 
is very simple to implement in existing MoM codes that are 
based on Galerkin’s method [61]. Here, we restrict the analysis 
to surface currents in free space. The corresponding stored 
energies for dielectrics and lossy media are more involved 
and still not well understood [38]. 

In antenna current optimization, the currents include both 
the sources and/or excitation coefficients. These are then used 
to analyze array antennas, array-pattern synthesis and small 
antennas, etc. Wheeler [75] considered an idealized current 
sheet to analyze array antennas. In array synthesis the optimal 
performance of the beamwidth and sidelobe level [53, 55, 69], 
for instance, is used as the optimization parameter to determine 
the array excitation. It is assumed that the excitations for 
different elements can be specified arbitrarily and that these 
excitations generate the desired current distribution as well as 
the radiated field. This procedure has been very successful 
in the areas of radar and communication. Another example, 
namely superdirectivity [39, 68] can occur if the excitation is 
optimized for maximal directivity [59]. Moreover, the direc- 
tivity is unbounded for finite apertures. These superdirective 
arrays are impractical, however, since, the magnitudes of the 
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excitations are large, implying high losses and strong reactive 
near fields [39, 68]. The superdirective solutions in the opti- 
mization problems are avoided by incorporating constraints on 
the losses and the reactive fields [31, 55]. 

The radiation properties of antennas are considered in 
antenna current optimization. It is assumed that the current 
distribution can be controlled in the antenna region, meaning 
that the amplitude and phase of the current density can be 
prescribed arbitrarily in this region. Optimization is used to 
synthesize current densities that are optimal with respect to 
antenna parameters such as the Q-factor, gain, directivity and 
efficiency. It should be noted that the current density is in 
general non-unique for optimal performance [31, 35]. 

In this paper, a review of antennas, stored energy and convex 
optimization for antenna current optimization is presented. In 
particular, convex quantities in antenna analysis and electro- 
magnetics and their relationship to optimization are discussed. 
MATLAB codes for maximization of the gain to Q-factor 
quotient, minimization of the Q-factor for superdirectivity and 
antennas with a prescribed radiated field are provided. The 
MATLAB codes can be copied from the pdf-file and are also 
available for download. The convex optimization problems are 
solved using CVX [24, 25] and standard MATLAB functions. 
The provided codes and data can be used to construct the 
results presented in the paper. 

The remaining part of this paper is organized as follows. 
Basic antenna parameters are reviewed in Sec. II. Optimal 
antenna design and antenna current optimization is discussed 
in Sec. III. Expressions for the stored energy, Q-factors and 
bandwidth are given in Sec. IV along with their corresponding 
matrix formulations in Sec. V. Convex optimization and con- 
vex quantities in electromagnetics are discussed in Sec. VI. In 
Sec. VII, antennas are analyzed using convex optimization. A 
dual formulation for the G/Q problem is given in Sec. VIII. 
Generalized eigenvalue problems and their relation to the 
stored energies are presented in Sec. X and some concluding 
remarks are presented in Sec. XI. Appendices containing a 
table of notation, discussion of stored energy, discussion of 
non-negative energy, a derivation of the dual problem and 
MoM data are in Apps A, B, C, D and E, respectively. 


II. ANTENNAS 


Antennas are ’the part of a transmitting or receiving system 
that is designed to radiate or receive electromagnetic waves’ 
according to the IEEE standard [46] (see Figs | and 2). A 
transmitting antenna must be matched to the feed structure 
such that the transmitted power is accepted by the antenna. 
The mismatch is quantified by the reflection coefficient. We 
introduce the antenna input impedance to separate the feed line 
from the antenna. In many cases we have a transmission line 
with real-valued characteristic impedance. This requires the 
antenna to be self-resonant, i.e., to have a negligible reactance 
and a resistance close to the characteristic impedance. The 
antenna input impedance, Zin, can be written 


2P} + 4jw(Wm — We) 
Fin? 





Zin = Rin + jXin = ’ (1) 





Fig. 1. Visualization of a radiating capacitive loaded dipole antenna. The 
orange shapes are the spherical capacitive caps of the dipole and the contour 
colors represent electric field strength, where red is high field strength and 
blue is low field strength. 


where we have also used the time average power and stored 
energy in the lumped circuit elements to express the input 
impedance [76], with the angular frequency w, dissipated 
power Pa, stored electric energy We, stored magnetic energy 
Wm, current Iin and imaginary unit j = /—1. 

There are many other parameters characterizing the perfor- 
mance of antennas, such as those presented in [2, 46, 74]. 
These are: 


bandwidth fə — fı and fractional bandwidth B = (fo — 
fı)/fo, where [f1, f2] is the frequency interval where 
the antenna performs according to the requirements and 
fo = (f2 + fi)/2 is the center frequency. The bandwidth 
requirements are usually formulated in terms of matching 
and radiation properties. 

directivity D(7) is the ratio of the radiation intensity in a 
direction 7 to the average radiation intensity [46]. The 
partial directivity denoted D(7,é) includes the depen- 
dence on the polarization ê. 

gain G7) is the ratio of the radiation intensity in the direction 
r to the average radiation intensity that would be obtained 
if the power accepted by the antenna were radiated [46]. 
The partial gain denoted G(7, ê) includes the dependence 
on the polarization é. 

efficiency ner, defined as the quotient between the radiated 
power and the accepted power, relates the gain and 
directivity G = neg D. 

radiation patterns are either specified with the magnitude of 
the electric field created by an antenna, |F'(7)|, or with 
the polarization, amplitude and phase F'(#) of the far 
field F. 

specific absorption rate (SAR) quantifies the amount of 
power absorbed per mass of tissue. 


The antenna parameter requirements are application-specific. 
For mobile phones, we strive for a large bandwidth, high 
efficiency, low directivity and low SAR in the considered com- 
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Pg Near field region ^> (radiated field) 
(reactive and radiated field) 








Fig. 2. Reactive and radiated fields from a current density J (r) in the region 
2 [37]. The reactive fields are concentrated in the near-field region around 
Q and vanish far from the source region 2, where E(r) ~ F(#)e—J*" /r. 


munication bands. The requirements for base station antennas 
are similar to those of mobile phones for bandwidth, but have 
usually higher efficiency and directivity. 

For small antennas, we often reformulate the fractional 
bandwidth in terms of the Q-factor [48, 74, 79]. The Q-factor 
is defined as the quotient between the time-average stored 
energy and dissipated energy 


2w max{W., Wm} 

Pa 
where we have also introduced the electric and magnetic 
Q-factors Qe = 2wW./Pa and Qm = 2wWm/Pa (see 
also Sec. IV). Radiation properties of antennas are equally 
important as matching, and are often used to characterize 
antennas (see Figs 1 and 2). The electromagnetic fields are 
also useful to determine the antenna Q-factor [12, 15, 29, 30, 
37, 58, 70, 79]. The radiated electromagnetic field is generated 
by oscillating currents on the antenna structure (see Fig. 2). 
From the radiation point of view, we can even consider antenna 
design as an art to produce the desired current distribution to 
achieve the radiation specifications. This can be thought of 
as simply modifying the antenna by shaping its structure and 
choosing its material properties to obtain the desired current 
distribution. 


Q= 





= max{Qe, Qm}, (2) 


III. OPTIMAL ANTENNA DESIGN AND CURRENT 
OPTIMIZATION 


Design requirements on antennas are often formulated in 
terms of combinations of antenna parameters such as those 
introduced in Sec. II. In addition to these parameters the 
antenna design is restricted by its size, weight and price where 
it is often desired to have small size, low weight and low 
cost. This often leads to contradictory goals as; for instance, 
electrically small antennas have narrow bandwidths and low 
directivity (D ~ 1.5) [12, 27, 32, 33, 37, 65, 71, 74, 77, 78]. 
Therefore, optimization is used to trade antenna performance 
versus size [44, 52, 62]. 

Antenna optimization is simply, optimizing the antenna 
structure with respect to the antenna performance in a given 


a) Maximal size of the antenna 


22K 


lx 
PO 


b) Antenna geometry with feed point 


[O O L 


c) Current distribution on the antenna 


[L L L 


d) Current distribution in the antenna region 








Fig. 3. Antenna and current optimization. In antenna optimization, we design 
antennas with optimal performance. In current optimization, we synthesize 
current densities with optimal performance. a) maximal antenna region. b) 
possible antenna design. c) current density on the antenna structure. d) possible 
current density in the antenna region. 


design space. Consider an antenna region defined by a planar 
rectangle, (2 = 924, with width ¢, and height @,, as depicted 
in Fig. 3a. The antenna optimization problem is to design 
an optimum antenna with respect to some parameters in 
the region (2 = Na by proper placement of metal (PEC) 
and feed. Fig. 3b depicts a center-fed PEC meander line 
antenna (one of many possible antenna designs) that fits within 
Na. The corresponding current distribution is depicted in 
Fig. 3c. The antenna current optimization problem is to find 
an optimal current distribution in the region 924, in this case 
for radiation in the broadside direction (see Fig. 3d). The 
obtained current distribution is not restricted to any specific 
feed-point or other constraints in the region (2,4. This implies 
that the current distribution from any antenna in the region 
Qa is a possible candidate for the current distribution in the 
design space. Consequently, the optimal current can be used 
to determine physical bounds (fundamental limitations) for 
antennas restricted to the design region 924. 


The previous example is sometimes encountered in practice 
but, antenna designers are often requested to design antennas 
in (small) parts of devices such as mobile phones, laptops and 
sensors. This case is illustrated in Fig. 4 where the structure 
N is divided into two regions; an antenna region Qa C 2 
and the remaining part Ra = Q \ Ra. Here, we refer to 
Na as the ground plane although it can in principle be any 
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Fig. 4. Device geometry with a region (2 with current density J, (left) with 
a device and (right) the geometry. We assume that the currents Ja can be 
controlled in the antenna region (24. The currents Jg in Ra = 2\ Qa are 
induced by the currents J, (see also [37]). 


type of region (metal or dielectric). Consider a typical device 
geometry to illustrate the approach (see Fig 4). The device 
structure is denoted by {2 and consists of an antenna region 2a 
and other components such as screen, battery and electronics. 
We assume that the antenna designer is allowed to specify the 
spatial distribution of the metal and dielectrics in the antenna 
region (2,4. The electromagnetic properties of the remaining 
region Ra = 2\ Qa are assumed to be fixed. For the antenna 
current optimization, we assume that the current density J A 
in (2, is controllable and that the current density Ja in Ra 
is induced by Ja. 

We can now formulate several optimization problems. The 
basic case with minimal Q-factor can be stated as: 


ee stored energy 
minimize = 2w———— (3) 

radiated power 
for lossless antennas. We minimize the Q-factor for an antenna 
by changing the material properties in the antenna region 2, 
for fixed material properties in Ng. For the approach in this 
paper, it is advantageous to rewrite the optimization problem as 
a constrained optimization problem. The minimal Q-factor (3) 
is then reformulated as minimization of the stored energy 

subject to a fixed radiated power P, = Pyo, i.e., 


minimize stored energy 


(4) 


subject to radiated power = Po. 


The two formulations (3) and (4) are equivalent but the latter 
formulation is more powerful as it is easily generalized by 
adding additional constraints. Alternatively, the optimization 
for the Q-factor can be formulated as maximization of the 
radiated power for a fixed stored energy, 


maximize radiated power 


subject to stored electric energy < Wo (5) 


stored magnetic energy < Wo, 


where W is a fixed number and the equality constraints are 
relaxed to inequalities. The relaxed problem contains stored 
energy equal to Wo as a special case, hence, the solution to (5) 
always gives a larger or equal radiated power than the problem 
with an equality constraint for the stored energy. Moreover, at 
least one of the inequality constraints in (5) will always be an 
equality (i.e., an active constraint) as otherwise the radiated 
power could be increased. 

We can easily generalize the optimization problem to many 
other relevant antenna cases. The quotient G/Q between the 
partial gain G = G(?f,é) and the Q-factor is investigated 


in [12, 32, 33]. The G/Q quotient gives a balance between 
a desired (high) gain and a low Q-factor. The G/Q problem 
can be written [31, 35] as 


minimize stored energy 


6 

subject to partial radiation intensity = Po. (6) 

Using that the partial radiation intensity is the squared mag- 

nitude of the far field [2, 5], we can rewrite the equality 

constraint in (6) into a linear equality constraint [31, 35]. This 
yields the optimization problem [31], which can be stated 


minimize stored energy 


i (7) 
subject to farfield = Fo. 


Antenna optimization problems can be solved by different 
methods that are suboptimal due to the numerical complexity 
of antenna problems. We can characterize the optimization 
approaches as local, global, model-based and combinations 
thereof. Local or gradient-based optimization is used to im- 
prove the design [16, 42, 43, 50]. This works very well if 
the initial design is close to the optimum otherwise there is a 
risk of getting trapped in a local suboptimal design. Global, 
stochastic, or metaheuristic optimization algorithms, such as 
Genetic Algorithms [44, 62], particle swarm [63], simulated 
annealing and Monte Carlo, are often used. These methods 
are very general and can be applied to almost any object 
functional. Model-based optimization and combinations of 
local and global algorithms can also be used [52]. For antenna 
current optimization, the problem is relaxed to determining the 
optimal current distribution instead of designing the antenna. 
Therefore, these problems can often be formulated as convex 
optimization problems, and can hence be solved efficiently [6, 
24, 31]. 

Below we illustrate the antenna current optimization and the 
associated physical bounds for the cases with Na = 9 and 
Na C N (see also Figs 3 and 4). The region (2 is a planar 
rectangle with side-length /x and ly (see Fig. 3a). The physical 
bounds are compared with data from classical dipoles, folded 
dipoles, loops and meanderline antennas in Sec. III-A and data 
from Genetic-Algorithm-optimized antennas in Sec. III-B. The 
results indicate that there are antennas that perform close to 
the physical bounds. This suggests that the antenna current 
optimization can be used to make an a priori estimate of the 
optimal antenna performance. 


A. Example: current optimization and physical bounds 


A planar rectangular structure is used to illustrate the 
problem of antenna current optimization for the G/Q bound 
in (6). The rectangles are infinitely thin and have the length 
lx and widths 4, = {0.5,0.1,0.01}4¢,. The quotient between 
the partial gain and the Q-factor G(Z, @)/Q is maximized for 
radiation in the normal direction of the plane (2-direction) 
and the polarization £. Fig. 5 depicts the upper bound on 
G/Q for ¢/X < 0.5 (half-a-wavelength). The bounds are 
identical to the forward scattering bound! [32, 33] for small 
structures [27]. The bound on G/Q improves with increasing 
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Fig. 5. Upper bounds on G(z,&)/Q for rectangular plates with height Zx 
and widths 2y = {0.5,0.1,0.01}2x, for 2x/A < 0.5, polarization ê = & 
and radiation in the * = 2 direction. G/Q from simulations of PEC strip 
dipole, capacitive dipole, meander and folded meander antennas are included 
for comparison with the physical bounds. The antenna feeds are indicated with 
a dot. The size of the dipole is ĉy = 0.01£x, the other antenna dimensions 
are ly = {0.5, 0.1}, with different line widths. 


antenna size l, and electrical size ¢,/A. This is a result 
of extending the degrees-of-freedom of the currents on the 
structure. 

The physical bounds are then compared with numerical 
results for self-resonant dipoles, folded dipoles, loops, me- 
anderline and folded meanderline antennas. The antennas are 
simulated in the commercial electromagnetic solver FEKO [1]. 
All of the antennas are matched to 50 Q input impedance and 
the antenna Q-factors are determined from (9). The simulated 
antennas have G/@ quotients close to the physical bounds 
(see also [3, 27, 32, 64] for additional comparisons). The strip 
dipole is resonant around y = 0.47 and has an optimal 
performance according to the G/Q metric. It can also be seen 
from the meanderline and folded meanderline antennas that 
the G/Q performance increases with antenna thickness. On 
the other hand, the resonance frequencies are shifted up as the 
effective length of the antenna is increased. 


B. Example: Genetic Algorithm and Current Optimization 


Antenna current optimization can be combined with global 
optimization algorithms. The former optimization is used to 
determine the physical bound on an antenna parameter, e.g., 
the Q-factor, directivity, radiation pattern, etc. The latter 
optimization is used to synthesize structures that perform opti- 
mally. Initial investigations of this automated optimal antenna 
design is considered in [13] and [14] for single- and multiband 
antennas, respectively. 

Here we maximize the partial-gain-@Q-factor quotient G/Q 
in (7) for electrical dimensions £/À < 0.5. The G(Z,«%)/Q 
quotient is considered for the 2-direction and #-polarization. 
The structures are considered infinitely thin perfect electrical 
conductors (PEC). They are restricted to rectangular regions 
in the xy-plane with the length £ = x and width 4, = £/2. 
The physical bounds on G(Z, @)/Q are computed with convex 
optimization (see Sec. VII) and depicted by solid lines in 
Fig. 6. The bounds are computed for an antenna restricted 


f / GHz, £ = 10cm 0/2 
0.9 1.2 1.5 


0.01 








Fig. 6. Solid lines—physical bounds on G(2,)/Q for antenna regions 
Qa restricted to rectangular regions. 6%, 10% and 25% of the region at 
the upper end in the x-direction is used for Qa, cf., Fig. 4. The situation 
with the entire region Qa = (2 used for optimization (100%) is included 
for comparison. Marks—G(Z,#)/Q values of structures optimized using 
a genetic algorithm (GA) [13, 14]. Insert—illustration of the considered 
situations. Blue and gold colored regions are the antenna 924 and ground 
plane g regions used in the convex and GA optimization (see also [37]). 


to 6%, 10% and 25% of the region at one end in the ¢,- 
direction (see insert in Fig. 6). Also, the case when the entire 
rectangular region (i.e, 100%) is used for optimization is 
included for comparison. The former three cases have been 
used in a Genetic Algorithm (GA) to synthesize the antennas. 
The G/Q quotients obtained by the GA-optimized structures 
are indicated as marks in Fig. 6. The presented results show 
that GA-synthesized antennas perform close to the physical 
bounds of the analyzed situations. 


IV. STORED ENERGY, Q-FACTOR AND BANDWIDTH 


The Q-factor is a measure of losses in a system and a high 
Q-factor describes a system with low losses. An oscillator with 
a high Q-factor will oscillate for a long time after the excitation 
has been removed. In antenna applications we want to dissipate 
power out from the antenna (see Fig. 2), hence a low Q-factor 
is desired. The Q-factor for an antenna tuned to resonance 
is defined as the ratio between the maximum of the stored 
electric (We) and magnetic (Wm) energies and the dissipated 
power (2). The electric and magnetic Q-factors correspond to 
the stored energy in the capacitors and inductors, respectively, 
normalized by the dissipated power in the resistors for lumped 
circuit networks. The time average stored energy in capacitors 
and inductors are 




















Civ |z/? LI |v 
M= ~ -E and Was — H -Er 
k V 
= | ME i a 
T po 

C 


respectively. Synthesis of lumped circuit networks leads to 
an alternative method to estimate the Q-factor from the input 
impedance of antennas [29]. 

The fractional bandwidth is inversely proportional to the 
Q-factor, i.e., a high Q-factor implies a narrow bandwidth. 
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Fig. 7. Magnitude of the reflection coefficient |I| for RLC circuits with 
resonance frequency wo and Q-factors Q = {6, 10, 30} [37]. The fractional 
bandwidths (8) for the Q = 6 case with threshold levels To = {1/V2, 1/3} 
are B ~ {0.33, 0.12}. 


The precise proportionality depends on the shape of the 
reflection coefficient. We can often quantify this shape with 
the distribution of resonances. The simplest case of a single 
resonance corresponds to series or parallel RLC circuits 


C oe 1, LR 
eG 


where the fractional bandwidth for single resonances is [79] 


gad i 2 
QA- Q 


and I denotes the threshold of the reflection coefficient. The 
reflection coefficients for single resonance RLC circuits with 
Q = {6, 10, 30} are depicted in Fig. 7. 

The estimate (8) is very accurate for Q >> 2 for the RLC 
circuit. The special case of the half-power bandwidth B ~ 
2/Q predicts an infinite bandwidth for Q = 1. This suggests 
that the Q-factor is most useful for Q >> 1 and in practice it 
is often sufficiently adequate for predicting the bandwidth if 
Q>5or Q > 10. Of course, the bandwidth can be increased 
by using matching networks [17]. 

Differentiation of the input impedance Z;n is a practical way 
to approximate the Q-factor for antennas [36, 49, 79] 


w| Zi, in Zl wR, 2 + wX in a Xin 2 
by _ VORP FOX Xn? gy 
oe 2Rin 2Rin 














for Ip = 1/V2 (8) 














where Zin,m denotes the input impedance tuned to resonance 
with a series capacitor or inductor. The formula (9) is exact 
for the series RLC single resonance circuits and often very 
accurate for antennas with Q >> 1 but can underestimate the 
Q-factor for lower values of Q, where multiple resonances 
are common [29, 36, 66]. For accurate estimates, (9) requires 
that the first order derivative |Z7, m| (linear term) dominates 
over the second- and higher-order derivatives. The relationship 
between the fractional bandwidth and Q-factor (8) for the RLC 
resonance circuit can also be used to define an equivalent Q- 


factor for a given threshold level I ie., 
2 Io 
Brn, a 


where Bp, denotes the fractional bandwidth for the threshold 
Io. 

In this paper, we estimate the Q-factor for antennas using 
the differentiated input impedance (9) and the Q-factor QB 
from Brune synthesized lumped circuit models [7, 29, 76]. 
The estimated Q-factors are used to compare the performance 
of antennas with the derived physical bounds from current 
optimization (see Figs 5 and 6). 

To analyze the radiation properties of antennas, we need to 
express the stored energy in terms of electromagnetic fields or 
current densities (see Fig. 2). The total time-harmonic energy 
is unbounded due to a large contribution from the radiated 
field far from the antenna [72, 73] for the corresponding time- 
domain case. This radiated field does not contribute to the 
stored energy of the antenna and is subtracted from the total 
energy [15, 21, 29, 58, 73, 79]. In this paper, we restrict the 
analysis to currents in free space (see also [38]). 

The integral expressions by Vandenbosch [70] represent the 
stored energy as quadratic forms in the current density (see 
also Geyi [21] for the case of electrically small antennas). 
The expressions are particularly useful as the radiated fields 
are generated by the current density on the antenna structure 
and hence directly applicable to the problem of current op- 
timization [31, 35]. The integral expressions are identical to 
subtraction of the energy density of the radiated far field for 
many cases [30] (see also App. B). 

The stored electric and magnetic energies are given by [30, 
70] 


Qn = (10) 


M cos(k 
We a Bf fvg (rı) )V2- J (r 2) ——— ( riz) 


Jd Arkri. 
3 . „~ sin(krj2) 
= (k Jı: J3 —Vı-JıV2° J3) — 5- dV1dV2 (11) 
and 
5 : cos(krj2) 
=X Jfk Iw = l ere 4rkriz 
22 
44 k 
= (k? Ji -J3 — Vi- Jı V2: 7) ee) ay, dV2, (12) 
T 


respectively, where Jn = J (rn) for n = 1,2, rig = |r1— rol, 
the asterisk * denotes the complex conjugate and we note that 
no /W = o/k. We also have the radiated power [22, 29, 70] 


Ra=B ff (sey) Te) 


22 


_ a : J(r1)V2 : J*(r2)) sin(k|r, = 


r2l) 


—_———- dV, dV2. 
4rk|rı — rol i 2? 


(13) 


For the radiation pattern and the directivity, we use the 
radiated far field [5, 60], F(#) = re!*" E(r) as r = |r| — oo, 
in the direction (see Fig. 2). The far field for the polarization 
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ê and direction 7 (with r - ê = 0) can be expressed as 


* A —jk Ax ikp. 
ê. F(*#) = Ee fe -T(rijel*F™ dy. aA 
T 
Q 
The partial radiation intensity is 
-p _ lê- FR? 
P(r,e) = (15) 
(7, é) Ona 
and the partial directivity and gain are 
An P(r, é€ An P(r, ê 
D(?,é) = ane and G(#,é) = ee (16) 


respectively, where P, + Po = Pa. In addition to the Q- 
factor (2), we consider the partial gain to Q-factor quotient 
G(r,é) dr P(7, è) _ nlé* - F(#)|? 
Q w max{W., Wm} wmo max{ We, Wm} 


that replaces the total radiated power in (2) with the radiation 
intensity. 








(17) 


V. MATRIX FORMULATION 


We consider a region {2 in which the current density J = 
J(r) is excited (see Figs 2, 3 and 4). This current density is 
expanded in local basis functions w,, as 


N 
I(r) = XO Inthn(r), (18) 

n=1 
where we introduce the N x 1 current matrix I with the 
elements J,, to simplify the notation. For simplicity, we also 
restrict the analysis to surface-current densities. The basis 
functions are assumed to be real-valued and divergence- 
conforming, with vanishing normal components at the bound- 
ary [61]. For simplicity, we use rectangular elements and basis 
functions with piecewise constant divergence (charge density) 
(see Fig. 8). Triangular elements with RWG or higher-order 
basis functions can also be used [61]. Moreover, we normalize 
the basis functions with their widths (cross-section for the vol- 
ume case) giving basis functions with the dimension length~ + 
(SI-unit m~t). The expansion coefficients are currents with 
the SI-unit ampere (A) and the impedance matrix (19) is in 
ohm (Q). It is easy to use dimensionless quantities by a scaling 
with the free space impedance no. 

A method of moments (MoM) type implementation using 
the Galerkin procedure is used to compute the energies given 
in (11) and (12). A standard MoM implementation of the EFIE 
using the Galerkin procedure obtains the impedance matrix 
Z=R+jX 


Zoom =m | | (ebm): balra) 
RQ 
+ zy Pm (r1)V2 -Yn (r2))G(r1 — r2)dS1 dS2, (19) 


e`jkr 


where the Green’s function [5] is G(r) = S~ and r = |r|. 
The expansion coefficients I are determined from ZI = V, 
where V is a column matrix with the excitation coeffi- 
cients [61]. 








Fig. 8. Illustration of discretization for a region (2 using rectangular mesh 
elements. The region is divided into the antenna region (24) and ground 
plane region, Qg. The amplitudes of six basis functions (18), three (green) 
in (2, and three (red) in Ng, are depicted. We let overlapping basis functions 
belong to the antenna part. 


Differentiating the MoM impedance matrix with respect to 
the wavenumber k gives 
kOZmn 
no Ok 





A 1 

= f | (ebm Yaa = Yr Praz bna 

QL 

+ (Kpm Una — Vi WmiV2° Wn2)ri2) G12 dS; dS2, 

(20) 

where Gi2 = G(r1 — r2), Yni = YnlT1), Pno = Vn(r2), 
n =1,...,N and rig = |rı — r2|. The MoM approximation 
of the stored energies (11) and (12) can be written as 








ly f/OX X 1-H 
We ~ gi (Z =) I= z! X.I (21) 
for the stored electric energy and 
ly /0X X 1 oy 
Wm ~ gl (+2 )r-51 Xml (22) 


for the stored magnetic energy, where the electric X, and 
magnetic Xm, reactance matrices are introduced and the super- 
script 4 denotes the Hermitian transpose. The expressions (21) 
and (22) are identical to the stored energy expression (for 
surface current densities and free space) introduced by Van- 
denbosch [70] and were already considered by Harrington 
and Mautz [40]. The total radiated power (13) for a lossless 
structure can be written as the quadratic form 


1 
P. x zI RI with R = Re{Z}. (23) 


We note that the computation of the reactance matrices and 
radiation matrix only require minor modifications of existing 
MoM codes. This makes it very simple to compute the stored 
energies and the additional computational cost is very low 
compared to the overall MoM implementation. Using the 
reactance matrices X. and Xm the EFIE impedance matrix 
is expressed as 


Z=R+j(Xm - Xe), (24) 
where we also note the relation 
IZ x 2P} + 4wj(Win — We) (25) 
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that resembles the energy identity for the input impedance (1). 
The Q-factor (2) for an antenna tuned to resonance can 
be expressed using the reactance matrices and the radiation 
resistance matrix 


2w max{W., Wm} 


q= P. + Pa 





_ max{IĦX.I, I"X,,1} wI" XI + [TAXI | 5i 
THR, + Re)I  2E(R,+Ro)I ’ pol 
where R = R,+Ro and Po = IĦRolI is the power dissipated 
due to ohmic losses. In this paper, we restrict the results to 
lossless structures so Ro = 0 and R = R, (see also [28, 34, 
38]). 
The far field (14) projected on ê is approximated by the 
N x 1 matrix FI ~ ê* - F(#) defined as 


= ETS ee 


Inserting (21), (22) and on in (14) we express the partial 
gain to Q-factor quotient (17) as 


G(7,é) 4n|FI|? 
Q mm max{ X.I, HXI} 
The electric and magnetic near fields [60] are approximated 
using the matrices Ne and Nm defined from 


Dana- Yo tom [5 Vi: 


= jk, (71)G (r = rı) dSı 








T Tı 





(r 1 ) ds 1: (27) 








(28) 


n(T1)VG(r — rı) 


(29) 


and 


N 
H(r) =NnIl= >> In fae) x ViG(r—r1)dS1, (30) 


respectively, where r ¢ 2. 

Embedded antennas (see Figs 4 and 8) are modeled with an 
antenna region where we can control the currents and a sur- 
rounding structure (ground plane) with induced currents [13, 
14, 31]. For simplicity, we restrict the discussion to induced 
currents on PEC ground planes. The induced currents depend 
linearly on the currents in the antenna region, and we use the 
EFIE (19) to determine the linear relationship 


Zaa Zac\ (la) _ (Va 

Zaca Zac) (Ie 0 
between the currents as discussed in [13, 14, 31]. The first 
row is unknown but the second row represents a constraint 


Zcala + Zecelge = CI = 0 (32) 


which can be added to the convex optimization problems in 
this paper. The decomposition of the basis functions into its 
antenna (Iq) and ground plane (Ig) parts is non-trivial as each 
basis function is supported on two elements (see Fig. 8). Here, 
we let the basis functions with support in both 24 and Ng 
belong to the antenna part I,. 

In the following, we assume that the numerical approxima- 
tion is sufficiently accurate so the approximate equal to (~) 
in (26) to (30) can be replaced with equalities. 


(31) 


not convex 





flax + By) 


Fig. 9. Convex and non-convex functions. Convex functions satisfy f (ax + 
By) < af(x) + f(y) for a+ 8 = 1 and, a,8 > 0, ie., the curve is 
below the straight line between two points (see [6] for details). Note that the 
non-convex function g is convex if the domain is restricted to the left or right 
of the local maximum in the middle. 


VI. CONVEX OPTIMIZATION AND CONVEX QUANTITIES IN 
ELECTROMAGNETICS 


Convex optimization problems are solved with efficient 
standard algorithms (see e.g., [6, 19, 24]). There is no prob- 
lem with getting trapped in a local minimum since a local 
minimum is also a global minimum [6] (see Fig. 9). A 
convex optimization problem is also associated with a dual 
problem. Dual problems are used to obtain posteriori error 
estimates. When an optimization problem is formulated as a 
convex optimization problem it is considered to be solved. 
There are of course difficult convex optimization problems and 
they can for example be ill-conditioned. Linear programming 
(LP), quadratic programing (QP) and quadratically-constrained 
quadratic programing (QCQP) are special cases of convex 
optimization. 

Convex functions f : RY — R satisfy [6] 


flax + By) < af(x) + Bfly) 


for all aœ, € R,a+f = 1, a,8 > 0 and x,y in the 
domain of definition of f. A simple interpretation is that 
the curve is below the straight line between two points for 
convex functions (see Fig. 9). Smooth convex functions have 
a positive semidefinite Hessian, i.e., the N x N matrix H with 


(33) 





elements H;; = soe. For functions of a single variable, 
the Hessian simplifies” to a non-negative second derivative 
af = f(x) > 0. A simple example is the second-order 
amal 


f(x) =ar? +br +c (34) 


that is convex if a > 0, as seen from the fact that f” (x) = 2a. 
A function g(a) is called concave if —g(a) is convex. The 
linear function f(x) = bx is both convex and concave. 

In this paper, we mainly use the following convex functions 


linear form f(x) = bx for 1 x N matrices b. 

quadratic form f(x) = x'Ax for symmetric positive 
semidefinite N x N matrices A > 0. 

norms f(x) = ||Ax|| 

max max{ f(x), fo(x)} for convex functions f(x), fo(x) 

logarithms — log(x). 
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We follow the convention in [6] and consider convex 
optimization problems of the form 


minimize f(x) 
subject to g(x) <0, t=1,...,m (35) 
Ax=b 


where the functions f(x) and g;(x) are convex and A is 
a matrix. In convex optimization, we can minimize convex 
quantities and maximize concave quantities. The linear (affine) 
quantities are both convex and concave so they can be either 
minimized or maximized. 

In order to study complex-valued quantities (e.g., electro- 
magnetic fields), we need to extend the definition of convexity 
to complex-valued functions. This can be achieved by con- 
sidering the real and imaginary parts as separate real-valued 
quantities. For our case, we note in particular that Re{-} and 
Im{-} are linear operators and that quadratic forms for positive 
semi-definite real-valued symmetric matrices A = AT are 
convex in the real and imaginary parts, i.e., 


z4Az = (x+jy)PA(x+jy)=x?Ax+y'Ay. (36) 


Convex optimization offers many possibilities to analyze radi- 
ating structures in terms of the current density. The expansion 
of the current densities in local basis functions (18) and the 
corresponding matrix approximations for the stored energy, 
radiated power and radiated fields are simple matrix operators 
on the current (see Sec. V). 

Examples of quantities commonly found in electromagnet- 
ics that are linear, quadratic, normed and logarithmic in the 
current matrix I defined in (18) are: 
linear: near fields N.I (29) and NmI (30), far field FI (27) 

and induced currents CI (32). 
quadratic: radiated power I"R., absorbed power, stored 
electric energy EIX. (21), stored magnetic energy 
ZI" XI (22), ohmic losses $1 Rol. 
norms: field strengths ||NI||2, far-field levels ||FI||2. 
max: stored energy for tuned antennas W = max{ We, Wm}. 
logarithmic: channel capacity. 





We can in general minimize convex quantities and hence 
convex optimization is very powerful to minimize (or restrict 
the amplitude of) power and energy quantities such as the 
stored energy, ohmic losses, radiated power, radiation intensity 
and side-lobe levels. This agrees with the goal of antenna 
design with the exceptions of radiated power. Consider for 
example the problem of minimization of the Q-factor in (5) 
where we have a finite stored energy (a convex constraint) 
but we maximize the radiated power. This is not a convex 
optimization problem as we should minimize convex quanti- 
ties. The corresponding minimization of the radiated power is 
convex and has the trivial solution 0 for I = 0. The same 
problem appears to apply to the gain Q-factor quotient (G/Q) 
in (6) and (28), where we minimize the stored energy for 
a fixed (partial) radiation intensity. This G/Q problem can 
however be reformulated to a fixed far field (7) that is linear 
and hence both convex and concave. In the following sections, 
we first illustrate the G/Q formulation and then generalize the 
formulation to superdirectivity and embedded antennas. 


VII. CONVEX OPTIMIZATION FOR ANTENNA ANALYSIS 


Optimization can be used to determine optimal currents 
and physical bounds for many relevant antenna problems [31, 
35]. Convex optimization offers great flexibility to analyze 
and formulate optimization problems [6, 31] and is directly 
applicable to G/Q in (28). Maximization of the partial gain 
to Q-factor quotient is analyzed in Sec. VII-A, applied to 
strip dipoles in Sec. VII-B and implemented using CVX in 
Sec. VII-C. Minimization of the Q-factor for superdirective 
antennas is considered in Sec. VII-D. Short dipoles and 
embedded antennas are analyzed in Secs VII-E and VII-F, 
respectively. Relaxation and a dual formulation is used refor- 
mulate the G'/@-problem in Secs VII-G and VIII, respectively. 


A. Partial gain to Q-factor quotient 
The partial gain to Q-factor quotient (28) in the used MoM 
approximation (18) is bounded by maximization of (28) over 
the current matrix, i.e., 
G(r, é) Z 4r|FI|? 
———— < max : 
Q  ~ I  max{I*X,I, EXI} 
Using the scale invariance of G&/Q with respect to I, namely 
that G/Q is invariant to the complex scaling of I > al, we 


can rewrite the maximization of G'/Q into minimization of the 
stored energy for a fixed partial radiation intensity 


max{T"X.I, HX mI} 
|FI|? =1, 





(37) 


minimize 
(38) 

subject to 
where the dimensionless normalization |FI|? = 1, or equiv- 
alently |FI| = 1, has been used. Moreover, the scaling 
invariance shows that we can consider an arbitrary phase 
FI = —j that removes the absolute value [31]. The particular 
choice used here is due to the presence of —j in (14) and 
produces real valued currents on planar structures for maximal 
radiation in the normal direction. In total, we get the convex 
optimization problem to minimize the stored energy for a fixed 
far-field in one direction and polarization [31] 


minimize max{I"X.I, I" X mI} 


; , (39) 
subject to FI = -j. 


Let I, denote a current matrix that solves (39). The minimum 
value of the stored energy in (39) is unique although the 
current vector I, is not necessarily unique. The optimum 
solution yields an upper bound on G/Q for the considered 
direction 7 and polarization ê, i.e., we have 
G(?f, è) 2 G(?, è) _ An|FI,|? 
Q E Q opt E no max{I} Xelo, HX mIo} l 
(40) 
The convex optimization problem (39) can be rewritten as 
follows. A normalized stored energy w = 4wW is introduced 
to obtain the equivalent (convex optimization) formulation 








minimize w 

subject to ha a He w, 
IXI < w, 
FI = -j. 


(41) 
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Fig. 10. A thin strip dipole with dimensions £x, ly divided into Nx = 16 
rectangular mesh elements. Two piecewise linear divergence conforming basis 
functions are depicted. 


The formulation (41) is referred to herein as the primal 
problem (P) (see App. D). An alternative optimization formu- 
lation is also to maximize the far field for a bounded stored 
energy [31]. 


B. Example: strip dipole 


We consider a planar rectangular structure to illustrate the 
antenna current optimization and physical bounds on G/Q 
in (40). The rectangle is infinitely thin and has length ¢ = £x 
and width 4, = 0.02¢, (see Fig. 10). The G/Q is maximized 
by (39) for radiation in the normal direction of the plane (2) 
and polarization 2. 

To maximize G/Q, we first compute the electric reactance 
matrix X. and magnetic reactance matrix Xm from (21) 
and (22). We can use local basis functions on triangular 
elements, rectangular elements or global basis functions; such 
as trigonometric functions. In this example, we start with a 
rather coarse discretization using Ny x Ny = 16 x 1 identical 
rectangular elements (see Fig. 10). The translational symmetry 
gives Toeplitz matrices 


X. = toeplitz(X.1) and Xm = toeplitz(Xm1) (42) 


where Xe1 denotes the first row of Xe and correspondingly for 
Xmi. The far-field matrix F is an imaginary valued constant 
column matrix. Altogether we have the MATLAB code below 
to compute the above quantities: 


oO 


Parameters and data for a 0.48\lambda strip dipole 


eta0 = 299792458 x* 4e-Txpi; % free space impedance 

kl = 0.48 x 2xpi; % wavenumber, 0.48lambda 

Nx = 16; % number of elements 

N = Nx-1; % number of unknowns 

dx = 1/Nx; % rectangle length 

dy = 0.02; % rectangle width 

Xell = le3*[1.14 -0.4485 -0.0926 -0.0153 -0.0059 
-0.0030 -0.0018 -0.0013 -0.0009 -0.0008 
-0.0007 -0.0006 -0.0005 -0.0005 -0.0004]; 

Xe = toeplitz(Xell); % E-energy 

Xm11 = 10*[1.8230 0.8708 0.2922 0.1664 
0.1060 0.0680 0.0411 0.0208 0.0050 
-0.0074 -0.0171 -0.0244 -0.0297 -0.0332 
=0.0351] 7 

Xm = toeplitz(Xm11); % M-energy 

Rr11 = 0.1*[7.0919 7.0668 6.9918 6.8680 6.6974 





6.4824 6.2264 5.9331 5.6067 5.2521 4.8744 
4.4788 4.0707 3.6558 3.2393]; 
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Rr 
F 


toeplitz (Rr11)+eye (N)»2e-5; 
eta0« (-lixkl) /4/pixones (1,N) *dx; 


2 
© 


far field 


for a strip dipole with length 4x = 0.48, or equivalently ké, = 
0.48 - 27 ~ 3. Here, we use a fixed numerical precision to 
simplify notation. Also, the radiation resistance matrix is made 
positive semidefinite by the addition of a small diagonal matrix 
(see App. C). More accurate values and refined discretizations 
are considered in App. E-A. 


Dd 
~~ 


C. CVX implementation 


There are several efficient implementations that solve con- 
vex optimization problems. Here we use CVX [24], which 
provides the MATLAB code 


2 
© 


CVX code for maximization of G/Q 
cvx_begin 
variable I(N) 
variable w; 
minimize w 
subject to 


= 
6 


current 
n. stored energy 


complex; 


oe 


quad_form(I,Xe) <= w; % n. stored E energy 
quad_form(I,Xm) <= w; % n. stored M energy 
FxI == -1li; % far-field 

cvx_end 


oe 


GoQ = 4x«pi/ (wxeta0) 
x linspace(0,1,N+2); 
plot (x,real([0; I/dy; 


bound on G/Q 

x coordinates 

I/dy; 0])) 

for the maximization of G/Q in (40) using (41). CVX solves 
the convex optimization problem iteratively (see the CVX 
manual [24] for details) and gives G/Q =~ 0.3. This is 
consistent with a half-wave dipole that is self-resonant at 
L = 0.48, and the forward scattering bound D/Q < 0.3 
in [26, 32]. The resulting radiation intensity (15), radiated 
power (23), directivity (16), stored electric energy (21), stored 
magnetic energy (22) and Q-factors (2) for the resulting 
current distribution, are computed as 


= 2 
6 


0]),x,imag([0; 


oe 


antenna parameters from the max. G/Q problem 


P = abs(Fx*I)*abs(F*I)/2/eta0; % radiation intensity 
Pr = real(I'«*«RrxI)/2; % radiated power 

D = 4*pixP/Pr % res. directivity 

We = real (I'«XexI) /4/kl1; % stored E energy 

Wm = real (I'«Xm*I)/4/k1; % stored M energy 

W = max (We, Wm) ; % stored energy 

Q = 2*klxW/Pr % Q 

Qe = 2«klxWe/Pr; % Q electric 

Qm = 2«*klxWm/Pr; % Q magnetic 


The normalized electric and magnetic stored energies are Qe ~ 
Qm * 5 and the directivity is D ~ 1.65, for the strip dipole 
data in Sec. VII-B. 

The current density distribution is depicted in Fig. 11. Here, 
we note that the current density J = Jy®& + Jyy is real- 
valued. This is because the special case with 7 = 2 gives an 
imaginary-valued far field vector F and, hence, a real-valued 
current matrix I as seen from (36). This a priori knowledge 
can be used in the CVX formulation above by declaring 


a 
© 


variable I(N); real valued current 


In this presentation we continue to use the complex-valued 
form to simplify the notation and avoid errors when we treat 
the general case with 7 F 2. 
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Fig. 11. The optimized current distribution on the strip dipole with length 
£ and width £/50 discretized with Nx = {16, 32} rectangles in the blue and 
red curves for the half wavelength case £/X = 0.48 (wavenumber k£ ~ 3). 
The radiation pattern, with D(2, &) ~ 1.64, is also depicted. 


We also note that it is preferable to use ||Xo 


I*X.I to replace the quadratic forms quad_form (I, Xe) 
with norms norm(sqrtXe*I) in CVX [24], where 
sqrtXe=sqrtm (Xe), which yields the modified MATLAB 
code 


2 
/ I|)? 


o 


% CVX code for maximization of G/Q 

sqrtXe sqrtm (Xe); 

sqrtxXm sqrtm (Xm); 

cvx_begin 

variable I(N) 

variable w; 

minimize w 

subject to 
norm (sqrtXexI) 
norm(sqrtXm*T) 


current 
sqrt stored energy 


complex; 


ole 


oe 


<= wW}; sqrt stored E energy 


sqrt stored M energy 


oe 


<= w}; 


Fel == -1li; > far-field 
cvx_end 
W = Ww; % n. stored energy 


oe 


GoQ = 4«pi/ (wxeta0) 
Pr = real(I'«*«RrxI)/2; 


bound on G/Q 
radiated power 


oe 


D = 2*pi/Pr/eta0; % directivity 
Q = w/Pr/2; sO 
x = linspace(0,1,N+2); % x coordinates 


plot (x,real([0; I/dy; 0]),x,imag([0; I/dy; 0])) 

where we have used the fact that the radiation intensity (15) 
is P = |FI|?/(2m) = 1/(2ņo) due to the normalization 
FI = -j of the far field in the optimization problem (41). 
The reformulation with norms improves the convergence, but 
requires pre-computation of the matrix square roots. We have 
observed that CVX works well for reasonably-sized problems 
and additionally solves the dual problem for improved per- 
formance [6] (see also Sec. VIII). Similar to Example VIU-B 
it is also important to make sure that the reactance matrices 
Xe and Xm are symmetric and positive semidefinite [35] (see 


App. C). 


D. Superdirective antennas 


Superdirective antennas have a higher directivity than a 
typical antenna of the same size [4, 39, 51, 57]. The directivity 


11 


given by (16) hints that the partial directivity is at least Do if 


_ Ane”. FF 


JP p 2 mle” FE 
2No P; l 


Do < D 
vo E no Do 


=> P, (43) 


This is added as the convex constraint $1"R,I < 27/(nDo) 
to the optimization problem (39) as follows 





minimize max{I"X,I,1'X,,1} 
subject to FI = -j (44) 
ERI < — 
no Do 


with the CVX code 


CVX code for minimization of Q for D\geq DO 
DO = 2; % directivity 
cvx_begin 

variable I(N) 
variable w; 
minimize w 
subject to 


oe 


current 
stored energy 


complex; 


ole 


quad_form(I,Xe) <= w; % stored E energy 

quad_form(I,Xm) <= w; % stored M energy 

imag(F*I) == -1; % far-field 

quad_form(I,Rr) <= 4*pi/DO/eta0;%radiated power 
cvx_end 


oe 


GoQ = 4x«pi/ (wxeta0) ; bound on G/Q 


Pr = quad_form(I,Rr)/2; % radiated power 

D = 2*xpi/Pr/eta0; % res. directivity 
Q = w/Pr/2; % res. Q 

x = linspace(0,1,N+2); % x coordinates 
plot (x, real ([0; I/dy; 0]),x,imag([0; I/dy; 0])) 


where we also note that the quadratic forms can be rewritten 
as norms for improved computational efficiency [24]. The 
resulting current is depicted in Fig. 12, where we observe 
the typical sub wavelength oscillatory current distribution for 
superdirective antennas [57]. The Q-factor is increased to 
Q = Qe ~ 160 for D = 2 in comparison with Q ~ 5 for the 
G/Q case (39) with D ~ 1.65. Moreover, the discretization 
Nx = 16 we used is not sufficient for an accurate description 
of the current. The case with N, = 32 is added and reduces 
the Q-factor to Q ~ 150. 


* J,(2) 








-0.5 





Fig. 12. Current distribution on the strip dipole with length £% = £ and width 
£/50 discretized with Nx = {16, 32} rectangles for the half wavelength case 
£/ = 0.48 (wavenumber k£ ~ 3). The radiation pattern, with D(2, &) ~ 2, 
is also depicted. 
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E. Short dipole 


Reducing the size of a dipole conserves the shape of the 
radiation pattern but adversely affects the Q-factor and thus 
the bandwidth. We consider a short dipole by increasing the 
wavelength to \ = 10¢. The MATLAB code is 


So 


Parameters and data fora 
= 299792458«4e-Txpi; 


.1\lambda strip dipole 
free space impedance 
wavenumber, 


oe 


oe 


Nx % number of elements 

N = Nx-1; % number of unknowns 

dx = 1/Nx; % rectangle length 

dy = 0.02; % rectangle width 

Xell = le3*[5.4722 -2.1527 -0.4441 -0.0729 
=0..0272 =0.0133 -0.0075 —0.0046 -0.0031 
=0.0022 -0:0016 -0:0012 =0:0009 -0:0007 
-0.0006]; 

Xe = toeplitz(Xell); % E-energy 

Xm11 = [3.8082 1.8348 0.6484 0.4050 0.2968 


0.2340 0.1926 0.1630 0.1407 0.1232 0.1091 
0.0975 0.0876 0.0792 0.0718]; 
Xm = toeplitz(Xm11); % M-energy 
Rril le-2«* [3.0819 3.0815 3.0800 3.0777 3.0743 
3.0701 3.0649 3.0587 3.0516 3.0436 3.0347 
3.0248 3.0140 3.0024 2.9898]; 
= toeplitz (Rr11)+eye (N) *3e-6; 
eta0« (-li*xkl) /4/pixones (1,N) *dx; 





Rr 
FS 


2 
© 


far field 


with the CVX code listed in 
Sec. VIL-C, gives G/Q 0.0028 and assuming a lossless 
structure Q ~ 544 and D = 1.5. In this case the stored electric 
energy Q = Qe ~ 544 (540 with Nx = 32) dominates over 
the stored magnetic energy Qm œ% 25, i.e., the short dipole 
is capacitive. This dipole illustrates the design difficulties of 
small antennas; reduced size also reduces the bandwidth. The 
current distribution of the short dipole is similar to that of a 
half-wave dipole antenna (see Fig. 11), since the dipole mode 
is relatively invariant under decreasing antenna length [35]. 


Solving the short dipole 


w 
xX 


w 
~ 


F Embedded antennas 





Fig. 13. A thin strip dipole with dimensions £x, y divided into Nx = 16 
rectangular mesh elements. The region is decomposed into the antenna region 
Qa in the center with 4 elements and the surrounding PEC ground plane 
region (2g at the edges with totally 12 elements. The corresponding basis 
functions are divided into 5 basis functions in the antenna region (2, and 10 
basis functions in Ng. Two basis functions are depicted, one in 24 (green) 
and one in Ng (red). 


Current optimization is easily generalized to the case of 
antennas embedded into a (PEC) ground plane as depicted 
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in Fig. 4. The current density is then decomposed into the 
controllable current I, and the induced current Ig. Similarly, 
the region is divided into (2,4, the antenna structure and Qc, 
the ground plane. In 24 we can fully control the currents I,, 
that in turn induce Ig in Qe [13, 31]. We use that Maxwell’s 
equations are linear implying that Ig depends linearly on 
Ia (32). 

The strip dipole geometry with length @ = ¢, and width 
fy is decomposed into the antenna region N4 in the center 
and ground plane regions 2g at the edges (see Fig. 13). 
This resembles a center-fed strip dipole with an extended feed 
region. The coupling matrix C in (32) is computed using the 
MATLAB code 


% antenna region in the center of a strip dipole 

Nf = 7; % start of antenna region 
indA = [Nf: (N-Nf+1)]; % antenna region indices 
indG = [1:(Nf-1) (N-Nf+2):N];% ground plane indices 
zm = Rr+lix (Xm-Xe); % EFIE impedance matrix 

Cm = Zm(indG,:); % induced current Cm*I=0 


The constraint CI = 0 (32) is added to the G/Q optimization 
problem (41) giving the convex optimization problem 


minimize w 

subject to Ix.I< w, 
IXI < w, (45) 
FI = -j, 
CI=0 


with the corresponding CVX code 


2 


% max. G/Q for an embedded antenna structure 
sqrtXe = sqrtm (Xe); 
sqrtXm = sqrtm (Xm); 


cvx_begin 
variable I(N) 
variable w; 
minimize w 
subject to 
norm(sqrtXex1) 
norm(sqrtXm*1) 


oe 


current 
sqrt stored energy 


complex; 


oe 


oe 


<= w}; 
<= w}; 


sqrt stored E energy 
sqrt stored M energy 


oe 


F*I == -1li; % far-field 

Cm*«I == 0; % induced currents 
cvx_end 
w = W*w; % n. stored energy 


oe 


GoQ = 4x«pi/ (wxeta0) bound on G/Q 


Pr = real(I'«*RrxI)/2; % radiated power 
D = 2*pi/Pr/eta0; % directivity 
Q = w/Pr/2; % 0 


Zz 


© 


x linspace (0,1,N+2); x coordinates 

plot (x,real([0; I/dy; 0]),x,imag([0; I/dy; 
0]),x(1+indA) , 0*x(indA),'d') 

The constraint CI = 0 can alternatively be used to eliminate 

Ic from the optimization problem [31]. 

The G/Q quotient for the short dipole case in Sec. VII-E is 
considered with the discretization N, = {16,32}. The center 
of the strip is used for the antenna (feed) region (2,4. The first 
case has 2 elements and the second case has 10 elements for 
Nx = 16 and twice as many for Nx = 32. This corresponds 
to widths of 24/16 = 0.125 and 100/16 = 0.6252 (see 
the circular marks in Fig. 14). The obtained G/Q values 
are {0.0022, 0.0027} with the corresponding Q-factors Q ~ 
{677,551} for the two cases. Increasing the discretization to 
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Fig. 14. The optimized current density distribution on the strip dipole with 
length £ and width €/50 discretized with Nx = {16,32} rectangles in the 
blue and red curves for the case £/X = 0.1 (wavenumber k£ œ~ 0.63). 
The antenna region consists of the 2 (blue curve) and 10 (red curve) center 
elements in the Nx = 16 case and twice as many in the Nx = 32 case. The 
width of the antenna region is also marked by circles. 


Nx = 256 reduces the Q-factors slightly to Q ~ {673,546}. 
The resulting current density is depicted in Fig. 14. Here, 
it is seen that the current density approaches the triangular- 
shaped current distribution on a short center-fed dipole as Qa 
decreases [2] (see also [31, 35]). The antenna region can thus 
be considered as the feed region. 


G. Relaxation of G/Q and Pareto fronts 


Before we introduce the dual problem for the G/Q opti- 
mization (41), we note that the solution of (41) produces a 
(unique) minimum value of w and a minimizing current I. 
The current I can be used to investigate if the constraints 
in (41) are equalities or inequalities. This answers if the 
antenna performance is constrained by the stored electric or 
magnetic energies. To further investigate the dependence of 
the stored electric or magnetic energies, we can consider 
simultaneous minimization of W, and Wm. This leads to 
multicriterion optimization and Pareto fronts [6]. The Pareto 
front is determined using the scalarization 


1 


1 
aW.+(1—a)Wm = —I” (aX. +(1-a)Xm)I = —I®X,I 
4w 4w 
(46) 
with 0 < a < 1 and the optimization problem 
4r|FIa|? 
maximize EL (47) 


This problem has a closed form solution as shown below 
in (53) and provides information about the tradeoff between 
the stored electric and magnetic energies. Small values of a 
emphasize the stored magnetic energy, whereas large values 
emphasize the stored electric energy. 

The multicriterion optimization (47) simplifies the optimiza- 
tion problem (37) considerably as the max operator in the 
denominator is removed. This can alternatively be interpreted 
by use of the inequality max{ A, B} > aA + (1 — a)B for 
0 <a <1 to replace the max operator, i.e., 


max{I"X.1, HXI} > H(aX. + (1 — a)Xm)I = IXI 
(48) 
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in the denominator of (37). Maximization with this relaxation 
increases the value and gives an upper bound on G/Q 


G 4n|FI,|? 


Q opt no Xa Ia 


for 0 < a < 1. As the bound (49) is valid for all a, we can 
consider the minimization problem 

G 
Q 
associated with the maximization problem (37). This is an ex- 


ample of duality and, as shown below, provides an alternative 
way to solve (37). 





< max (49) 


a 


San 4r|FI,|? 
min he ee a 
T 0<a<1 nol#X ola 


(50) 





ma: 
Ia 
opt 


VIII. DUAL PROBLEM FORMULATION FOR G/Q 


When solving the optimization problems in Sec. VII-C, it 
is observed that CVX states that the dual problem is solved 
for improved efficiency. Duality is a powerful principle in 
optimization. Dual problems can be used to construct efficient 
algorithms, to estimate errors and to provide insight into the 
optimization problem [6]. 

The upper bound on G/Q in (37) is reformulated to the 
convex optimization problem (41). The primal problem (P) 
in (41) is associated with a dual optimization problem (D) (see 
App. D). The dual function d(a) (see the general definitions 
in (94) and (102)) is defined as the minimum value of the 
following optimization problem 


(aX. + (1 -— a)Xm)Ia 
subject to FI, = —j, 


minimize 


(51) 


where 0 < œ < 1 is a Lagrange multiplier. The explicit 
solution is given by (see (115) and (116)) 








4) = rex. 4 : a)Xm) F" ao a 
and 
: —1pH 
E o = SRF Sane 202) 
The dual optimization problem (D) is stated as 
maximize d(a) (54) 


subject to 0<a< 1, 


which can be solved efficiently by a line-search, such as the 
Newton’s method, bisection method, golden-section search and 
parabolic interpolation [19]. The dual formulation also yields 
useful interpretations and bounds in terms of the dual current 
Ia and the associated electric Qea, magnetic Qma and total 
Q-factors Qa = max{Qea, Qma } in (2), as well as the partial 
gain Ga. Here, weak duality (97) (dual cost < primal cost) 
implies that 


= 4T AQca + (1 ai a)Qma 

no Ga 
max{I!?X.I,, IH 
[FHT |? 





d(a) 


Xmla} = 4T Qa 


, (55) 
no Ga 
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which should be compared to the inequality (48). Use that I, 
is a suboptimal solution to (41) and rewrite in G/Q to get the 
final bounds 


Ga G Ga 
< < ; 
Qa Q opt aQea + (1 = a)Qma 


which is compared with the relaxation (50). 

Strong duality (the duality gap is zero) can be shown for 
the optimization problem (41) by using Slater’s constraint 
qualification [6] (see also App. D). Here, the strong duality 
implies that 





(56) 


š a G 
min = 
0<a<1 AQea + (1 m a)Qma Q opt 


It is noted that the dual formulation (54) normally provides a 
much more efficient way of maximizing the partial gain to Q- 
factor quotient defined in (37) than by direct use of the primal 
problem formulation (P) given by (41). 

The derivative of the dual function is given by 


odd dl 


~da dad 


where it is seen that d(a) increases (decreases) for capacitive 
(inductive) cases. The opposite is true for Go/(aQe + (1 — 
a)Qm), Le., the G/Q quotient decreases (increases) for ca- 
pacitive (inductive) cases. The second derivative of the dual 
function 





(57) 


d (Xe — Xm)Ia = -IEXIa, (58) 


„_ d id’? 

da? d 
can be used to solve (54) with Newton’s method based on the 
update 


=2 





xx x1, (59) 


Ant1 = An — d' (am) /d" (an). 


Newton’s method is very efficient for cases with the optimal 
value a, in the inner region 0 < 6 < ao < 1— ô such 
that d’(a,) = 0. Newton’s method can be combined with the 
bisection method or golden section search for cases with a, 
approaching 0 or 1. 

Small electric dipole type of antennas are often capacitive 
and have a dominant electric stored energy (We > Wm). This 
implies that d(a) is increasing (d’ > 0) and the maximal value 
of (54) is obtained for a, œ~ 1. The resulting Xa = Xe for 
a = 1 in (52) is often singular and hence difficult to invert. 
These cases with small electric dipole type can be solved with 
an initial a close to 1, e.g., a = 1—6 with ô = 1074. A 
simple algorithm would be to evaluate Ga /Qa and check if the 
duality gap is below the desired threshold. If not update 6 —> 
2*'6§ with the — sign (+ sign) for the capacitive (inductive) 
case, cf., the bisection method. 


(60) 





A. Numerical example: strip dipole 


The upper bound on G/Q for the strip dipole in Sec. VII-B 
is determined using the dual formulation (54) with the inequal- 
ity (56). The optimization (54) is solved using the MATLAB 
function fminbnd with the reformulation of max. d(a) as 
min. —d(a) together with the explicit solution (52). This 
yields the MATLAB code 
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Fig. 15. Gain Q-factor quotient, G/Q, from the dual optimization prob- 
lem (54) with the inequalities (56) for a strip dipole with size £ x €/50, 


/X = 0.47, 7 = Z and ê = ĉ. 








% maximization G/Q using a dual formulation 


da = @(a)-1/real (Fx ((axXet+(l-a) *Xm) \F'));% -d(a) 
[a,d] = fminbnd(da,0,1); % min -d(a) 
GoQ = -4xpi/eta0/d 


with the result ag ~ 1 and G/Q|opt ~ 0.3 for Nx = {16, 32} 
and £ = 0.47. The fminbnd function uses a combination 
of golden section-search and parabolic interpolation [19] and 
minimizes the functional until the error in œ is below a certain 
threshold. 

The inequality (56) can be used to obtain error estimates for 
G/Q|opt (for the used MoM approximation). The parameter a 
is used to determine the current I,, and the associated stored 
energies and Q-factors, i.e., 


% antenna parameters from the dual problem 
Xa = axXet+(1l-a) «Xm; 


J = Xa\F'; 

d = 1/real (F*J) % dual value 

Ta = -lixdsd; % current 

P = abs((Fx*Ia)*(Fx*la))*4*pi/eta0;% n. rad. int. 
QoGe = real (Ia'*Xexla) /P; % Qe/G 

QoGm = real(Ia'*Xm«Ia)/P; % Qm/G 

GoQa = 1/max (QoGe, QoGm) % G/Q 

GoQd = 1/ (axQoGe+ (1-a) *QoGm) % dual G/Q 

GoQd-GoQa % duality gap in G/Q 


yields the duality gap in G/Q|opt on the order of 1077 for the 
strip dipole. Note that this is an estimate of the error in the 
maximization of G/Q for the used numerical approximation, 
i.e., the MoM with Nx = {16,32}. It is essential to investigate 
the convergence of the MoM approximation by refinement of 
the discretization, i.e., to increase Ny and Ny. 

The inequality (56) for the N, = 32 case is depicted in 
Fig. 15. The function values used by fminbnd are depicted 
by the dots. Here, it is seen that the evaluation for a ~ 0.4 
already gives the bound 0.29 < G/Qlopt < 0.31. The gap 
is larger for smaller values of a and becomes negligible as 
a — 1. This is due to the weighting of the stored magnetic and 
electric energies in (51), which emphasize either the magnetic 
or the electric energy. 

The physical bound on G/Q for the strip dipole is depicted 
in Fig. 16 for £ < 0.5X. The strip dipole is divided into 
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Fig. 16. Gain Q-factor quotient, G/Q, from the optimization problem (39) 
for a strip dipole with sides £ x €/50, £/A < 0.5, # = Z, ê = & and 


Nx = {50,100}. The forward scattering bound [32, 33] on G/Q is also 
included. 





Fig. 17. A planar rectangular region with dimensions Zx x ly divided 
into Nx x Ny = 12 x 6 rectangular mesh elements. Four piecewise linear 
divergence conforming basis functions are depicted. 


Nx = {50,100} and Ny = {1,2} rectangular elements and 
the resulting Xe, Xm, R matrices are described in App. E-B. 
The dual formulation (54) is used to maximize G/Q. The 
result is also compared with the forward scattering bound on 
D/Q from the polarizability” and the generalized absorption 
efficiency 1/2 [27, 32, 33]. The differences between the 
results obtained with different discretizations and methods are 
negligible. 


B. Numerical example: planar rectangle 


Although the strip dipole geometry in Secs VII-B 
and VIII-A is very good for illustrating the optimization 
concepts it has a trivial polarization dependence and cannot 
radiate a magnetic dipole pattern efficiently (negligible loop 
currents). We consider a planar rectangle to obtain polariza- 
tion dependence and loop currents. We place the rectangle 
in the xy-plane and let the side lengths be &, = @ and 
ty = €/2 (see Fig. 17). To start, we consider an equidistant 
discretization using Nx = 2, Ny = 64 and hence a total of 
Nx(Ny — 1) + Ny(Nx — 1) = 4000 expansion coefficients 
(see App. E-E). This is a significant increase in optimization 


*http://www.mathworks.com/matlabcentral/fileexchange/26806-antennaq 
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Fig. 18. Gain Q-factor quotient, G/Q, from the optimization problem (39) 


for a rectangular plate with side lengths £ and £/2, wavelength A > 2¢, 
r= {%, 9, 2} and ê = {#, 9, (@+jy)/V2}. The G/Q is normalized with 
k3a3, where a = V5 /4 is the radius of the smallest circumscribing sphere. 
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Fig. 19. Resulting Q from the optimization problem (39) for a rectangular 
plate with side lengths £ and ¢/2, electrical size €/X < 0.5, directions * = 
{%,y, 2} and polarizations ê = {&, 9, (& + jĝ)/ V2} (see Fig. 18). 


variables compared to the strip-dipole case and it is also 
observed in the increased computational time to solve the 
optimization problems. 

The G/Q quotient is maximized for combinations of radi- 
ated fields in the 7 = {%, ĝ, Z}-directions and polarizations 
ê = {@, 9, (2+jĝ)/ V2}. The maximal gain Q-factor quotient, 
G/Q, normalized with k3a? is depicted in Fig. 18. The 
result is also compared with the forward scattering bound on 
D/Q from the polarizability? and the generalized absorption 
efficiency 7 = 1/2 [27, 32, 33]. The resulting Q-factor is 
depicted in Fig. 19 for £ < A/2, where we see that Q is 
lowest for the r = y direction and ê 2 polarization. 
The optimization (39) is solved using CVX [24] and using 
the dual formulation (54) with the fminbnd function in the 
MATLAB code and Newton iterations (60). The final results 
are indistinguishable but the Newton iteration is faster for 
larger problems. Note that several solvers can be used in 
CVX for improved performance (see [24] for details). There 
are also many quadratically-constrained quadratic program 
(QCQP) solvers with better performance. 


3http://www.mathworks.com/matlabcentral/fileexchange/26806-antennaq 
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Fig. 20. G/Q determined from the dual problem (54) and bound (56) for a 
rectangular plate with side lengths ¿x = 2y = 4/10, direction 7 = 2 and 
polarization ê = & (see Fig. 19). The maximal value is G/Q|opt ~ 0.0123 
giving Q = 125 and D & 1.53. 


The dual problem (54) is illustrated in Fig. 20 for the 
rectangular patch in Fig. 19 with £ = 0.1, and radiation in the 
fr = 2Z-direction for the ê = #-polarization. The four curves 
Ga/(aQea T (1 = a)Qma), Ga/ max{Qea, Quma}, Ga/Qea 
and Ga/Qma are depicted for 0 < a < 1. The stored electric 
energy dominates until a 1 and the resulting radiation 
pattern is similar to that of an electric dipole. We note that 
Ga/(AQea + (1 — a)Qma) decreases towards its minimum 
at a ~ 1, and contrarily, Gy/max{Qea,Qma} increases 
towards its maximum at a 1. The dual problem (54) is 
solved using Newton iterations (60) starting from ag = 0.5. 
The evaluation points are marked with circles in Fig. 20. The 
first iteration gives a; > 1 and then we set a; = 0.99 and 
combine the Newton and bisection methods. The convergence 
is very fast and the optimal value G/Q ~ 0.0123 is obtained 
after 3 iterations. The resulting current distribution gives 
Q 125 and D ~ 1.53. The corresponding results for 
Ny = 2Ny = 32 are G/Q ~ 0.0121, Q ~ 126 and D © 1.53. 
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Fig. 21. Resulting current density and charge density for the rectangular plate 
in Fig. (20). The current density is -directed and strongest at the edges. The 
charge density distribution is close to the charge density on a PEC plate in a 
static electric field [35]. 


The resulting current distribution and charge density p = 
ZV - J are depicted in Fig. 21 for the case N, = 32 and 
Ny = 16. The current density is aligned with the longest edges 
+a-directions) and concentrated close to the edges. The di- 
rection of the current density is however counterintuitive. The 
current is -directed in the edge elements but —a-directed in 
some neighboring elements. This is a small antenna structure 
£ = /10 that is dominated by the stored electric energy (see 
Fig. 20). The stored electric energy (11) can be approximated 





—~ 
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by the electrostatic energy [21, i‘ in the limit kf > 0, i.e., 


W. ~ ) as, dS2. 
deo dla ee 


Here, it is seen that the stored electric energy is determined 
by the charge density for small antennas (see Fig. 21). This 
charge density distribution is similar to the induced charge 
density on a PEC rectangle in an electrostatic field [35]. The 
corresponding current density is non-unique as V - (J + V x 
J.) = V -J for any Je. The term V x Je contributes to the 
magnetic energy and current densities of the form V x Je can 
be added without affecting max{W., Wm} as long as Wm < 
We (see also Sec. IX and App. C). 





(61) 
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Fig. 22. Ga/Qa from (56) for the plate in Fig. 19, 24x = 2ly = 0.1A, 


rf = y and ê = @. The resulting current distribution is depicted in Fig. 23. 
The radiation pattern is depicted for three values of a, see also Fig. 24. The 
maximal value is G/Q|opt ~ 0.0259 giving Q ~ 102 and D ~ 2.66. 


The corresponding case with radiation in the y-direction for 
the -polarization is depicted in Fig. 22. The stored energy 
is dominantly electric for low values of a, but it changes to 
dominantly magnetic when a œ~ 0.67 or above. This value 
of a also gives the maximum of Go/max{Qea;Qma} for 
the considered I,, and the minimum value of Ga/(aQea + 
(1 — a)Qma). The Newton iteration (60) converges 
as a œ {0.5, 0.73536, 0.67677, 0.66629, 0.66602, 0.66602} 
with the corresponding dual gap in G/Q approximately 
10~{?:2,3.48,16} | The optimal value is G/Q œ~ 0.0259 that 
results in Q ~ 102 and D ~ 2.66. The resulting current 
density is depicted in Fig. 23. The real part of the current 
density is an -directed current radiating as an <&-directed 
electric dipole mode. The imaginary part is a loop-type current 
density that radiates as a Z-directed magnetic dipole (see also 
Fig. 24). 


IX. MINIMUM Q FOR PRESCRIBED RADIATED FIELDS 


Maximization of G/Q aims to realize a low Q-factor and 
a large gain. The gain is related to the directivity by the 
efficiency G = nef D and the maximum directivity is in the 
range of 1.5 to 3 for small antennas. Hence, it is mainly the 
Q-factor that changes for small antennas (see Fig. 19). The 
Q-factor also increases rapidly if the antenna is excited for 
superdirectivity, as is seen to be the case in Sec. VII-D. 
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Fig. 23. Resulting current density for the G/Q problem in Fig. 22. Real and 
imaginary parts to the left and right, respectively. The real part is dominated by 
an £-directed current radiating as an -directed electric dipole. The imaginary 
part is a loop type current that radiates as a Z-directed magnetic dipole (see 
also Fig. 24). 


The obtained current distribution from the maximal G/Q 
problem can be used to compute a resulting Q-factor. This Q 
gives the lower bound on Q for small lossless antennas with 
dipole type radiation patterns, i.e., antennas with G = D = 
1.5 or G = D = 3. The directivity increases often with the 
electrical size of antennas, e.g., half-a-wavelength dipoles have 
D ~ 1.64 that is larger than D = 1.5 for the Hertzian dipole. 
Although, the resulting Q-factor from the G/Q problem is still 
a good estimate for bounds on Q, there is no guarantee that 
it is the lower bound on Q. 





Fig. 24. Illustration of spherical modes. Z-directed electric dipole (left) and 
Huygens source composed of an &-directed electric dipole and a Z-directed 
magnetic dipole (right). 


The G/Q problem can be reformulated to that of mini- 
mizing Q for a projection of the radiated field on the desired 
field [31] (see also [31] for other possibilities). Small antennas 
radiate as electric and magnetic dipoles and the radiation 
pattern of larger antennas can be described in spherical 
modes [5]. The optimization problem is identical to (41) with 
the change of F to the regular spherical modes expanded in 
basis functions given in (18) (see [31]). Here, we use the 
MATLAB function 


Fm for projection of spherical modes 

order of the mode, 1 for dipoles 
Fourier component (azimuthal), 0,1,..,L 
1 for TE and 2 for TM 

0 for even and 1 for odd 
phmodematrix(k,bas,meshp,[l mt s]); 


AP oP oO 


oe 


mo ct 3 mæ 


m 


and then either use CVX or the dual formulation to solve the 
convex optimization problem. 

Consider the planar rectangle with side lengths 4% = 24, = 
0.1, (see Fig. 17). The minimum Q-factor for radiation of an 
x-directed electric dipole mode gives Q ~ 120 with D ~ 1.5 
for the N, = 2, N, = 64 case. This can be compared with 
the G(Z, é)/Q case in Fig. 20 that has Q ~ 125 and D & 
1.53. The quotient G'/Q is approximately the same for the 
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two cases but Q and G = D are slightly lower for the case 
with a desired dipole mode. The corresponding case with the 
combined electric and magnetic dipole modes in Fig. 24 gives 
Q ~ 102 and D ~ 2.65. This is similar to the G(g, &)/Q case 
and the current density resembles the distribution in Fig. 23. 

The projection on the spherical modes can be interpreted 
as a minimization of the Q-factor where the radiated power 
is replaced with the radiated power in the considered mode. 
Consider a factorization of R, as R, = FĦF,, where F, is 
the far-field. The decomposition R, = FEF, is not unique 
and can e.g., be computed from a Cholesky decomposition of 
R, or a mode expansion. Here R, is first transformed to a 
positive semidefinite matrix (see Sec. C). The radiated power 
is rewritten as 


1 1 1 i= 
P, = ZIRI = S"F8F,1 = S|FsI? = 550 [Fs nI)? 
2 2 S 5| | 2 2| , | 


(62) 
where F, ,, denotes the n™ row of F,. The decomposition FI 
can be interpreted as a mode expansion of the radiated field. 
The lower bound of the Q-factor is the minimum of 


max{I"X,1, I XmI} 2 max{I"X, I, I" X mI} 
Ena |Foynll? [Fs,noll? 





(63) 


where Fs no is the far field of the desired radiation pattern. 
This optimization problem is mathematically identical to the 
G/Q problem in (37) if only one mode is considered or if R, 
is a rank-1 matrix. This problem can be solved with the convex 
optimization formulation in (41). Note that this is similar to 
the dual problem of (51) rewritten as the quotient 


IXI 


(64) 
which is a Rayleigh quotient with the rank 1 matrix FĦF in 
the denominator. 


X. EIGENVALUES 


We observe that the Q-factor (26) resembles a Rayleigh 
quotient that is efficiently analyzed using generalized eigen- 
values. However, the maximum of the stored energies in (26) 
is difficult to handle and has to be removed by explicitly 
assuming that either of the stored energies is larger. The 
G/Q quotient also has a closed form solution under similar 
assumptions [35]. 

Current optimization for the Q-factor (26) differs from the 
G/Q case (28) by the use of the radiated power instead of 
the radiation intensity. Although this difference appears to 
be negligible, minimization of the antenna Q is much more 
involved than maximization of G/Q. This is mainly due to 
the possibility to reformulate the partial radiation intensity 
|FI|? in the G/Q problem (28) as the field FI in (39) and 
hence obtain a convex optimization problem. That is, we can 
replace maximization of the radiation intensity (power) with 
maximization of the field strength (see also (63)). 
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A. Optimization for antenna Q 

Instead of convex optimization, we can use the fact that 
the Q-factor (26) resembles a generalized eigenvalue problem 
except for the max operator. A possible approach is to 
relax (26) using convex combinations of W. and Wm, Le., 


max{T#X.1, HXI} 








Q= ERI = max{Qe, Qm} 
(aX. —a)Xm)I 
> aQe+(1-a)Qn = ER 65) 


for 0 <a < 1, cf., (48). The lower bound on the Q-factor, 
Qib, is hence formulated as a minimization problem for the 
right-hand side of (65), i.e., a Rayleigh quotient that can be 
solved as a generalized eigenvalue problem 


(aXe + (1 = a)Xm)Io,n _ QanRlon (66) 


with the eigenvalues Qa, = AQea + (1 — @)Qma ordered 
ascendingly. Let I, denote the eigenvector associated with the 
smallest generalized eigenvalue (eigenmode) Qa,1 in (66) and 
Qea and Qma the corresponding electric and magnetic Q- 
factors. This gives the following estimate 


AQca + (1 = a)Qma < Qib < max{Qea, Quma} 


for the lower bound Qip. 

The solution of (66) is depicted in Fig. 25 for a planar 
rectangle with side lengths 4% = 244 = 0.1. The stored 
electric energy Qe dominates for a < 0.8 and the stored 
magnetic energy Qm dominates for a > 0.8. We note that 
the convex combination aQe + (1 — a@)Qm increases up to 
its maximum Qa 102 at a ~ 0.8. The corresponding 
Qa = max{Qea,; Qma} is close to its minimum Qa ~ 123 
for a range of 0.1 < a < 0.8. There is hence a minimum gap 
of approximately 21 between the minimal eigenvalue, Qa 1, 
in (66) and the realized Q-factor, i.e., we have the estimate 


(68) 


(67) 


Dd 
we 


102 < Oi, < 123 


for the lower bound Qi» of the considered region. 

One problem with the minimization of the Q-factor us- 
ing (65) and (66) is the lack of control of the radiated 
field. This is illustrated by the radiation patterns in Fig. 25. 
The electric dipole pattern dominates for the lower values of 
a, where the stored energy is electric. The pattern changes 
abruptly from an electric dipole pattern to a magnetic dipole 
type pattern around a œ~ 0.8, i.e., where the stored energy 
changes from electric to magnetic. 

The bound (67) can be combined with the convex opti- 
mization problems for maximization of G/Q in (39) and the 
minimization of the Q-factor for a desired radiated field in 
Sec. IX. The resulting Q-factor from the maximization of 
G(y,xz)/Q in Fig. 22 gives a resulting Q ~ 102 and as 
Qi» > 102 according to (68) Qi, ~ 102. 

The eigenmodes of (66) are orthogonal 


D mXalomn =T mRIan =0 (69) 


for m #n and Xa > 0 and R > 0. The modes form a basis 
if R is positive definite, R > 0. The resistance is however in 
general only positive semidefinite R = O (see also App. C). 
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Fig. 25. Q-factor from the optimization problem (65) for the rectangular 
patch in Fig. 19. The parameters in (67) are depicted. 


A self-resonant antenna has an equal amount of stored 
electric and magnetic energy We = Wm. This special case 
simplifies the Q-factor in (2) for a lossless antenna with 
a = 0.5 in (65), ie., 


Q _ w(W. + Wm) -_ Qe + Qm — 
= P, = 5 = 


I" (Xs+ Xm) 


I 
2ĦRI ` (70) 





2 


This is considered in [23, 40] and solved as the generalized 
eigenvalue problem 


(Xe + Xn)In = 20: ,RI,. (71) 


bm 
Here, it is essential to observe that the solution of (71) is in 
general not self-resonant, i.e., We + Wm. Moreover there is 
no simple relation between Qi and the fractional bandwidth 
for the untuned case. 


B. Characteristic modes 


Generalized eigenvalues are used to define characteristic 
modes of metallic structures [8, 9, 11, 20, 41, 54]. The EFIE 
impedance matrix (19) is used to formulate the generalized 
eigenvalue problem 


ZI =(R+jX)I=(1+jA)RI or XI=ARI. (72) 


The eigenvalues A are real valued and the eigenvalues with 
the smallest magnitude are most significant. The eigenvalue 
problem (72) is associated with the Rayleigh quotient 


o I¥XI _ xe, —xX,)I 
~ JARI THRI 

The Rayleigh quotient (73) resembles (70) but with the 
difference between the stored energies instead of the sum. 


The characteristic modes strive for a low reactive power 
(resonance) instead of a low stored energy. 





A (73) 


C. Reduction of the number of degrees of freedom 


Eigendecomposition of the X., Xm, R matrices can be used 
to reduce the number of unknowns in optimization problems. 
Consider one of the optimization problems in this paper, e.g., 
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(Xm—Xe)In = ARI,  (Xm+Xe)In = ARI, 


TITID) PDI I AAAS 
aa vy 


ne vy 
rere eee CSCS Gh RR ELE 





Fig. 26. Eigenmodes associated to the three smallest magnitudes of 
generalized eigenvalues of (Xm F Xe)In = AnI for a planar rectangle with 
side lengths Z and £/2 and wavelength A = 102. The rectangle is discretized 
with Nx = 2Ny = 32 equidistant elements. 





the G/Q problem (41). Assume that the resulting current 
distribution I has the Q-factor Q < Qo, i.e., 


I?x.1 Xl 

TR = % H 

IFRI IFRI 
This implies that it is interesting to consider the subspace of 
current matrices that satisfy these inequalities. Subtracting and 
adding the inequalities suggests the eigendecomposition 


hg 0 ee a, a 
IRI 

or equivalently to determine the eigenspace associated with 
the generalized eigenvalues A < Qo to (Xe F Xm)I = VRI. 
There is however no requirement that the solution to the 
optimization is an eigenmode. Assume for simplicity that the 
optimal current is the sum of two eigenmodes, i.e., I = I; +I2, 
with corresponding eigenvalues A, and As. Then the orthog- 
onality (69) implies 


and (74) 





< Qo. 





< 2Qo (75) 





(L F l2)" (Xe F Xm)(L + Le) 
(L + In)FR(I, + I2) 
_ TE(X. F Xm) + E(X. F Xm) 
~ ERT, + ERI 

















(76) 





and hence that a high Q-factor for mode I> does not imply a 
high Q-factor for I, F Iz because the denominator consists of 
the sum of the dissipated powers of the modes. 

The reactance matrices X. and Xm have very few or no 
negligible eigenvalues, i.e., they have full rank. The radiation 
resistance matrix R, has many small eigenvalues that can be 
discarded to reduce the number of unknowns in the optimiza- 
tion problem (degrees of freedom) (see Sec. C). Numerical 
tests indicate that it is more efficient to use the eigenspace 
induced by the generalized eigenvalues from (66) or (72), i.e., 





(Xm F Xe)In = A4 nin R (17) 











with the smallest magnitude |A+,n|. The smallest 45 eigen- 
values for the planar rectangle with 4% = 2¢, and 4% = 
{0.1, 0.25, 0.5}A are depicted in Fig. 27. There are potentially 
N = 4000 eigenvalues but it is only approximately 20 that are 
reliable due to the spectrum of R for l% = 0.1A (see Fig. 28). 
The eigenmodes for (71) and the characteristic modes (72) 
are similar. The first three eigenmodes are depicted in Fig. 26. 
Their radiation patterns are similar to the patterns of « and 
y-directed electric dipoles and a Z-directed magnetic dipole, 
respectively. The eigenvalues decrease as the electrical size 
increases (see the £x = {0.25,0.5}A cases). 




















1 5 10 15 20 2 30 35 40 


Fig. 27. The smallest 45 eigenvalues |A+,,,| of (77) for the planar rectangle 
with £ = lx = 24y and £ = {0.1,0.25,0.5}A. The rectangle is discretized 
with Nx = 2Ny = 64 equidistant elements (see Fig. 26). 





Consider the eigenvalue decomposition (77) and order the 
eigenvalues An in order of ascending magnitude. Divide the 
eigenvalues such that A,, < 6 for 1 < n < N; and A, > ô for 
n > Nı where ô is the chosen threshold level for the negligible 
eigenvalues. Let the columns of U consists of the eigenmodes 
I, for n = 1,...,N1 < N normalized as I,/,/ITRI,, 
where it is used that the eigenmodes I, are real valued. This 
decomposition reduces the number of unknowns. 


Ix UI (78) 


that gives the approximation 
1 1 eyoars ~z l eyo z 
W. ~x —IPX.I x —FUTX.UI = —I"X.ÍI (79) 
4w 4w 4w 
and similarly for Xm, R and F, i.e., 


Xm = UTX,,U,R = UTRU, and, F=FU. (80) 


Approximating the optimization problem (39) leads to the 
minimization problem 
minimize max{I"X,1,I"X,, I} 


Be (81) 
subject to FI = -j. 


and similarly for the other optimization problems in Sec. VII 
and VIII. This reduces the number of unknowns from N to 
Ni. 

Maximization of G/Q using (81) with the N, = 20 << N = 
4000 smallest eigenmodes (77) makes negligible differences 
in (39). It is even sufficient to use the three smallest modes 
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N, = 3 for relatively high accuracy. The reduction of the 
number of unknowns can be very efficient for the solution 
of complex optimization problems. The reduction can also 
provide physical insight from the interpretation of the char- 
acteristic modes [8, 9, 11, 20, 41, 54]. The computational 
advantage for the G/Q type optimization problem (39) is 
however limited as (39) can be solved with a few Newton 
steps (60); moreover the generalized eigenvalue decomposition 
can have higher computational cost. 


XI. DISCUSSION AND CONCLUSIONS 


A tutorial description of antenna current optimization has 
been presented. The presentation is intended to illustrate differ- 
ent possibilities with the approach. The included examples and 
data are chosen to illustrate the theory and be simple enough 
to stimulate investigations using MATLAB and CVX . It is also 
straightforward to convert the codes to other languages. 

Antenna current optimization can be used for many common 
antenna geometries. In this tutorial, we have focused on the 
case with antennas occupying the entire region, ie., Qa = 2 
(see Fig. 3). The case with a PEC ground plane is also 
discussed (see Fig. 4 and [13, 14, 31]). Generalization to 
antennas embedded in lossy media is considered in [28] and 
antennas above ground planes in [67]. Geometries filled with 
arbitrary inhomogeneous materials can also be analyzed using 
optimization of the equivalent electric and magnetic surface 
currents [5] for some cases [47]. 

There are many possible formulations for the antenna cur- 
rent optimization problem. This offers a large flexibility and 
possibilities to model many relevant antenna cases. The simple 
case with maximal G/Q leads to minimization of the stored 
energy for a fixed radiated field in one direction (41). The 
generalization to antennas with directivity D > Do is obtained 
by addition of a constraint to the total radiated power (44). The 
stored energy can also be minimized for a desired radiated 
field or by projection of the radiated field on the desired 
far field [31]. The case with antennas embedded in a lossy 
background medium is very different since there is no far 
field in the lossy case. It is however simple, instead, to include 
constraints on the near field [28]. It is also possible to impose 
constraints on the sidelobe level or radiation pattern in some 
directions, cf., the cases in array synthesis [53, 69]. 

Validation of the results against simulations and/or measure- 
ments is very important. The bounds on G/Q are compared 
with classical antennas in Fig. 5 and GA optimized antennas 
in Fig. 6 (see also [3, 13, 29, 32, 45, 64]). It is essential 
to compute the stored energy and Q-factors for the antennas 
accurately in these comparisons. The Qz; formula (9) is very 
useful for single resonance cases but it can underestimate the 
Q-factor for cases with multiple resonance [36, 49, 66, 79]. 
The stored energy in circuit models synthesized from the input 
impedance offers an alternative approximation of the stored 
energy [29]. The circuit models are synthesized using Brune 
synthesis [7] technique, which requires an analytic model (PR 
function) of the input impedance from zero frequency and up 
to the frequency of interest [29]. 

The computed current densities can be used for physical 
understanding. For instance, the case with the half-wavelength 
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strip-dipole in Fig. 11 is recognized as the classical cosinus- 
shaped current distribution. This shape is also close to op- 
timal for longer wavelengths [35]. The oscillatory current 
distribution of the superdirective dipole in Fig. 12 resembles 
the case with superdirective arrays. The current distributions 
for more complex structures are more difficult to visualize. 
Typical dipole and loop currents are seen on planar rectangles 
in [31]. Here, it is important to understand that the value of 
the objective functional (e.g., G/Q) is unique but there are in 
general many current distributions that gives this value. The 
same holds for the derived quantities such as; the resulting Q- 
factor and, directivity calculated from currents that minimize 
G/Q. 

The accuracy of the convex optimization solution is easily 
verified using the dual formulation (56); hence, it is not a 
major problem. The underlying accuracy of the MoM type 
discretization of the problem is however essential for the 
reliability of the computed results. Here, as for all MoM 
solutions, it is important to investigate the convergence of 
the discretization to see how the results depend on mesh 
refinement. Moreover, if it is known a priori that a specific 
mesh is sufficient to model all antennas, then the same mesh 
can be used for current optimization. 


The accuracy of the expressions for the stored energies (11) 
and (12) are also essential for antenna current optimization. 
It is known that (11) and (12) equal the sum of the stored 
energy defined by subtraction of the energy in the far field 
and a coordinate-dependent term [29, 30]. The coordinate 
dependence vanishes for small structures and also for struc- 
tures with a symmetric radiation pattern [30, 79]. The stored 
energies (11) and (12) also reduce to the classical stored 
energies in the static limit. However, the stored energies (11) 
and (12) can produce negative values for electrically large 
structures [35]. This questions the validity of (11) and (12) for 
larger structures. The expressions have been validated against 
the Qz formula [79] and circuit models for several antennas 
in [13, 14, 29]. The values agree for cases with large Q-factors 
but can disagree as Q approaches unity [29]. This coincides 
with the region where Q is a useful concept and can be used 
as an estimate for the fractional bandwidth. In this tutorial, 
we have restricted the size of the structures to approximately 
half-a-wavelength (ka ~ 1.5 to 2). This is much larger than 
the classical definitions of small antennas ka < 0.5 or ka < 1. 
For the planar rectangle it can lead to low Q-factors and hence 
questionable results when comparing the antenna performance, 
e.g., Q = 1 corresponds to an infinite bandwidth using (8). 
There is still no consensus of the stored energies for larger 
structures and for inhomogeneous materials, so much research 
remains in these areas. 
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APPENDIX A 
NOTATION 


Scalars are denoted with an italic font (f, F), vectors (in 
IR?) with a boldface italic font (f, F) and matrices with a 
boldface roman font (f, F). We consider time harmonic fields 
in free space with the time convention e/“*. 


APPENDIX B 
STORED ENERGY 


The stored energy expression (11) is motivated by the 
identity [30] 


FOP ay 


r2 


Him | BoP- 


ro +00 4 
Ir|<ro 


= Wet Woo = 7. mf fv: J {V2° 
22 


-JıV2- J3 


apee Jı- J3 


-JıVo2: a 


J: costae) 


tkry9 


o 


— (k? Jı- J} dVı dV2 


72 by (kr) dV, dV2, (82) 


where Ji(K) = (sin(x) — Kcos(K))/&? is a spherical Bessel 
function and F is the far-field (see Fig. 2). The identity (82) 
is valid for arbitrary current densities with support in a 
bounded region (2 radiating in free space, (see Fig. 2). The 
derivation of (82) is solely based on integral identities for 
the free space Green’s function and vector analysis [30]. The 
integral in the left-hand side is the difference between the 
electric energy density and the energy density of the far- 
field term [18, 23, 79]. The first integral in the right-hand 
side is coordinate independent and identical to the stored 
electric energy W, in (11) proposed by Vandenbosch [70]. The 
second term W,.9 contains the coordinate dependent factor 
r? —r3 = (r1—12)-(r1 +72) and Weco has the coordinate 
dependence 


Wea Wig = + I d-*#|F(#)|? dS» (83) 


|*|= 


for a shift of the coordinate system r — d-+r (see [30]), where 
the integration is over the unit sphere. The expression (12) 
for the stored magnetic energy is motivated by the analogous 
identity 


_ Ho 2 |F(*)|? 
lim > J |H(r)|* — om dV = Wm + Weco. (84) 
Ir|<ro 


Note that the identities (82) and (84) are valid for current 
densities with arbitrary frequency dependence and that they 
differ from the expressions in [23] (see also [10)]). 
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Speed of light, co = 1//€0k0 
impedance of free space, no = \/ Ho/€o 
permeability of free space, po = no /co 
permittivity of free space, co = 1/(noco) 
electric field 

magnetic field 

current density 


Jn = J (rn) for n= 1,2 
charge density, p = zV -J 
far field 


input impedance 

input resistance, Rin = Re Zin 

input reactance, Xin = Im Zin 

stored electric energy 

stored magnetic energy 

dissipated power 

radiated power 

ohmic losses 

Q-factor (2) 

electric Q-factor (2) 

magnetic Q-factor (2) 

Q from Z/, (9) 

Q from the fractional bandwidth and Tọ 

Q from Brune circuit [29] 

reflection coefficient, see Fig. 7 

threshold level for the reflection coefficient (see Fig. 7) 
directivity, also partial directivity D(7, ê) 

gain, also partial gain G(#, ê) 

position vector in R3 (see Fig. 2) 

magnitude of r, i.e., r = |r| (see Fig. 2) 
distance |r1 — r2| 

(unit) direction vector, i.e., P = r/r (see Fig. 2) 
(unit) polarization vector (see Fig. 2) 

source region (see Fig. 2) 

antenna region, (24 C §2 (see Fig. 4) 

ground plane region (see Fig. 4) 

side length of a rectangle, also lx, ly (see Fig. 3) 
frequency 

angular frequency w = 2r f 

wavenumber k = w/co, kno = wHo, k/no = weo 
wavelength A = co/f 

basis function (18) 

current matrix 

impedance matrix (19) 

resistance matrix, R = ReZ 

reactance matrix, X = Im Z 

electric reactance matrix (21) 

magnetic reactance matrix (22) 

far-field matrix (27) 

near-field matrices (29) and (30) 

induced currents matrix (32) 

current matrix in the solution of dual problems 
Q-factor for the current Ig, 

electric Q-factor for the current Ig 

magnetic Q-factor for the current I, 

convex combination Qa = @Qea + (1 — a)Qma 
gain for the current Ia 

free space Green’s function, G = e~J*I"| /(47|r]) 
imaginary unit, j? = —1 

complex conjugate, (a + jb)* = a — jb 
transpose 

Hermitian transpose 

positive definite, IZAI > 0 for all I 40 
positive semidefinite, FAI > 0 for all I 

unit vector, || = 1 

Lagrange multiplier, also v for matrices 

nabla operator 

volume element 

surface element 
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APPENDIX C 
NON-NEGATIVE STORED ENERGY 


The integral expressions for the stored energies are not 
positive semidefinite for all structures [35]. In [30], this is 
interpreted as an uncertainty of the stored energy due to the 
subtraction of the radiated power in the interior of the struc- 
ture. The convex optimization approach in this paper relies on 
having positive semidefinite quadratic forms. The expressions 
are observed to be positive semidefinite for sufficiently small 
structures but can be negative when the size is of the order of 
half-a-wavelength [35] (see also Figs 29 and 30). In practice 
there might be some small negative eigenvalues for smaller 
structures due to the finite numerical precision in the MoM 
approximation and the relatively large subspace with small 
eigenvalues. Note that the stored electric energy at statics has 
an infinite dimensional null space consisting of all solenoidal 
current densities, e.g., of the form V x A for some vector 
field A. The resistance matrix has a null space containing 
non-radiating sources [5], i.e., current densities of the form 


1 
jwHo 
for vector fields f = f(r) with compact support. 





J = (kf-VxVxf) (85) 
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Fig. 28. Eigenvalues An of Xe, Xm, R for a planar rectangle with side 
lengths 4% = 2y = 0.1\ divided into 64 x 32 elements (4000 basis 
functions). 


The eigenvalues Ag, Am,n and Arn of Xe, Xm and R, 
respectively, for a planar rectangle with side lengths / and 
£/2 and wavelength A = 10@ are depicted in Fig. 28. The 
rectangle is divided into 64 x 32 identical elements giving 
64 x 31 +63 x 32 = 4000 basis functions. The definite sign of 
the eigenvalues Am,n > 0 shows that Xm is positive definite 
Xm > 0. The electric reactance matrix is also positive definite 
Xe > 0 with approximately half of the eigenvalues Aen 
much larger than the remaining ones. The small eigenvalues 
belong to divergence free eigenmodes and they approach 0 
as £/\ — 0. The resistance matrix R should be positive 
semidefinite but the finite numerical accuracy of the evaluation 
of (19) makes R indefinite. In Fig. 29, it is seen that R 
has a few (~ 20) dominant large eigenvalues and lots of 
small eigenvalues. The small eigenvalues are of the order 
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1071? smaller than the dominant eigenvalues. These small 
eigenvalues are very sensitive to the numerical evaluation of 
the impedance matrix (19) and similar to the eigenvalues of a 
random matrix, cf., the MATLAB plot 


A=rand (1000); 
semilogy (abs (eig(AtA"'))) 


These random errors results in approximately 2000 small 
negative eigenvalues. The negative eigenvalues are marked 
with circles. 


103 | 


1072 J 


1077 | 


10712 L 











O 500 1000 1500 2000 2500 3000 3500 4000 


Fig. 29. Eigenvalues An of Xe, Xm, R for a planar rectangle with side 
lengths £ and £/2 divided into 64 x 32 elements (4000 basis functions) for 
the wavelength A = 2¢. The negative eigenvalues are marked by circles, i.e., 
An < 0. The eigenmode (current) to the negative eigenvalue Ae,4000 is also 
depicted, cf., with the loop current in [35]. 


The corresponding case with the wavelength A = 0.5¢ 
is depicted in Fig. 29. The positive eigenvalues Ay», > 0 
show that Xm is positive definite Xm > 0. Xe has one 
negative eigenvalue A. 4900 showing that Xe is indefinite. 
The corresponding eigenmode (eigenvector) is an equiphase 
loop current as depicted in the inset, see also the explicit 
construction in [35]. The resistance matrix R is indefinite due 
to the used numerical accuracy similar to the A = 10€ case 
in Fig. 28. The matrix R has a few more dominant large 
eigenvalues compared to the A = 10@ case as the number of 
radiating modes increases with the electrical size 4/A. 

The negative eigenvalue for Xe vanishes for longer wave- 
lengths and X. = 0 for electrically small structures. Fig. 30 il- 
lustrates maximal size of a planar rectangle such that X. = 0. 
The rectangle side lengths ¢, and 4, are normalized with 
the wavelength. The longest side of the rectangle is divided 
into 32 equidistant regions and the highest frequency with 
Xe > 0 is determined using Cholesky factorization as depicted 
by the blue curve marked with circles. The corresponding 
value with an indefinite X, is illustrated by the red curve 
marked with circles. The results are similar to discretization 
using 64 regions. The region below the curves is the region 
with positive semidefinite X.. The corresponding Q-factors 
determined from the forward scattering bound [32, 33] on 
D/Q assuming D = 1.5 are depicted with the contours for 
Q = {1,2,5,10,20,100}. The classical regions for small 
antennas ka < {0.5, 1} are shown with the dashed blue quarter 
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Fig. 30. Illustration of the maximal size of a planar rectangle with side 
lengths x and £y such that Xe is positive semidefinite Xe = 0. The blue 
(red) curve with circular marks illustrate the largest (smallest) case with 
Xe positive semidefinite (indefinite). The contours illustrate the region with 
Q > {1, 2,5, 10, 20, 100}, where the Q-factor is estimated from the forward 
scattering bound [32, 33] assuming an electric dipole pattern. The classical 
regions for small antennas ka < {0.5,1} are also depicted with the blue 
dashed quarter circles. 


circles. The region where Xe is indefinite corresponds to Q- 
factors below 2 and hence values where Q loses its meaning 
and there is in practice no restriction on the bandwidth (8). 

In this paper, we consider the stored energy as zero if the 
integral expressions are negative [31]. This is performed by 
an eigenvalue decomposition of the reactance matrices Xe and 
Xm and the resistance matrix R, e.g., 


X. = UA,U?, (86) 
where A, is a diagonal matrix containing the eigenvalues 
Aen. Negative eigenvalues are replaced by 0, i.e., Aen > 
max{ Aen, 0}, giving the electric reactance matrix 


X. > U max{ Ae, 0}UT (87) 


and similarly for Xm and R. 

Although this approach eliminates the problems with indefi- 
nite matrices, it is not entirely satisfactory. One minor problem 
is that the reactance Xm — Xe is changed. This is however 
easily solved by addition of the same quantity to both Xe 
and Xm. In this paper, we restrict the size of the antenna 
structure to approximately half-a-wavelength to mitigate the 
problem with negative stored energies. This also coincides 
with the typical range where the antenna performance is 
restricted by the Q-factor (bandwidth) (see Fig. 30). The 
energy expressions produce reliable results for some simple 
antennas for substantially larger structures [29], but much 
research remains to be done before we can draw any definite 
conclusions. 
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Diagonalization of the reactance matrix can be used to sepa- 
rate X into two positive semidefinite matrices X = X,—X_, 
where X, = 0 and X_ > O. The simplest case is to 
diagonalize X, i.e., 


X = UAUT = UA, UT —- UA_UT =X, — X_, (88) 


where A+ > 0. Note that the decomposition is non-unique and 
that any positive semidefinite matrix can be added to X+ and 
X_. The decomposition (88) resembles the decomposition of 
the reactance matrix into the electric and magnetic reactance 
matrices (24). 





APPENDIX D 
DUALITY IN CONVEX OPTIMIZATION 


A brief overview of duality in convex optimization, and its 
application to the minimization of the maximum of quadratic 
forms is given below (see also [6, 56]). 


A. Primal and dual problems 


Consider the primal optimization problem (P), which can 
be stated as 


minimize f(x) 
subject to g(x) <0, (89) 
Ax = b, 


where x € R”, f(x) € R is a convex function, g(x) € R” 
a vector of convex functions, A € R?*” and b € R’. Let X 
denote the affine space X = {x € R” |Ax = b}. 


Define the Lagrangian function 
L(x, v) = f(x) +" g(x), (90) 


where v € R” is a vector of Lagrange multipliers. It is readily 
seen that for any x € X 


and hence the primal optimization problem (P) in (89) is 
equivalent to the minimax problem 


f(z) 
+00 


if g(x) < 0, 


otherwise, ih) 


L — 
mar (x, v) 


(92) 


minimize max L(x, v). 
xEX v>0 
The dual optimization problem (D) is defined by interchanging 
the order of the minimization and maximization, and is hence 
defined by the problem 


ae d(v), (93) 
where d(v) is the dual function 
d(v) = min L(x, v), (94) 
and where v > 0. 
It is readily seen that for any x € X and vp > 0 
in L(x,v) < L(x, vp) < D(x 
min L(x,0) < L(x) < mar (x, v) (95) 
and hence that 
max min L(x, v) < min max L(x, v), (96) 


v>0 xEX xEX v>0 
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or 
< min 


(97) 
xe X,g(x)<0 


max d(v) F(x), 


v>0 
which is a statement of weak duality. 

Slater’s constraint qualification [6] states that if there exists 
an x € X such that g;(x) < 0 for i =1,...,m (x is strictly 
feasible) in (89), then strong duality holds for the convex 
optimization problem, i.e., 


maxd(v)= min x). 98 
ea) an te) (98) 
B. Minimizing the maximum of quadratic forms 
Consider the convex optimization problem 
aceon. ce T m 
minimize max{x A;x}j) 
{ U Jier (99) 


subject to Ax = b, 


where x € R”, A; € R”*” are symmetric positive semidef- 
inite matrices, A € R?*” and b € R?. This problem is 
equivalent to the following primal formulation (P) 





minimize y 
subject to xTA;x < Y, (100) 
Ax = b, 


which is an optimization problem over (x, y) € R”+t, Let X 
denote the affine space X = {x € R”|Ax = b}. The dual 
function d(v) is given by 


d(v)= mi s(x” AGx — 
(v) e a (xTA;x — y)} 
= : . xl A. 
= agin, {ud = 2 vi) + 3 vixTA;x}, (101) 
or 
$ T ; 
min x v;A; |x if y, = 1, 
oe A a 
—oo otherwise, 


where v € R”™ is a vector of real valued and non-negative 
Lagrange multipliers, v; > 0. 

The linearly constrained quadratic minimization problem 
above has an explicit solution with 


m =1 m =1 
x= (£ vas) ATIA (>. was) AT] b. 
i=1 i=1 
(103) 
The dual optimization problem (D) in (93) becomes 
m —1 E 
maximize bT | A 2 nas) AT b 
i=1 
m (104) 
subject to 5 vi =l, 
i=1 
Vi 2 0, 


which consists of an m — 1 dimensional search over explicit 
solutions to the quadratic minimization problems (based on 
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convex combinations of the matrices A;) as defined in (102). 
The dual problem (D) in (104) can sometimes be computa- 
tionally advantageous to solve in comparison to the primal 
problem (P) in (100), in particular when m is small and the 
system dimension n is large. 

It is readily seen that the convex optimization problem (P) in 
(100) satisfies Slater’s constraint qualification [6] by choosing 
an arbitrary x € X and a y € R such that 


xTA;x—y <0, (105) 


for alli = 1...,m, ie., an (x, y) exists that is strictly feasible. 
Hence, strong duality holds for this convex optimization 
problem. 


C. Linearly constrained quadratic optimization in complex 
variables 

Let I = I +jI, € C”*1 where IL, I; € R”*!. The complex 
gradient with respect to I* is defined by 





o 1/02 o 
= j ; 1 
oF 2 (sr tide) (108) 
It is noted that the condition oe = 0 is equivalent to the 


Cauchy-Riemann equations when f is holomorphic, and to the 














condition or = oF = 0 when f is real valued. The following 
differentiation rules are also readily verified 
ð = 
ge FI=0, 
IEF" = F", (107) 
(2) == 
IT IFRI = RI, 


where F € C!*” and Re C”*”, 
Consider now the linearly constrained quadratic optimiza- 
tion problem 
I"RI 
subject to FI=4g, 


minimize 


(108) 
where g € C is a constant. The corresponding Lagrange 
function is given by 

L(I, v) = IFRI + Re {v* (FI — g)}, (109) 
or 


L(I, v) = IRI + n Re {FI — g} + 4 Im {FI — g}, (110) 


where v = vr +jv; is the complex-valued Lagrange multiplier. 
Now, since 


ô Re{v*(FI— g)} 








Or* 
ð fv*(FI—g)+v("F"—9*)| vom 

-2f 5 =F , (11) 

the condition for optimality is 

o v 

L(I, v)=RI+-F"=0 112 
Lv) = RI+ F" =0, (112) 

and hence i 
I= -3R F", (113) 
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The multiplicator v is found from the constraint requirement Xel1 =1e3* [1.57397 -0.57065 -0.15929 -0.02964 
FI = g. Hence, -0.01124 -0.00552 -0.00314 -0.00197 -0.00133 
v -1nH -0.00094 -0.00070 -0.00055 -0.00044 -0.00036 
— JFR F =g. (114) -0.00031 -0.00026 -0.00023 -0.00021 -0.00019 
-0.00017 =0 -00016 =0:00015: -0.00014 -0.00013 
The optimal solution becomes -0.00013 -0.00012 -0.00011 -0.00011 -0.00010 
-0.00009 =0: 0000977; 
g -1H Xe = toeplitz (Xe11); % E-energy 
I= RFE Be (115) xm11 = [6.77879 3.65774 1.50119 0.93418 
0.66881 0.50954 0.40106 0.32110 0.25892 
and the corresponding minimum value is 0.20865 0.16685 0.13135 0.10073 0.07402 
0.05053 0.02977 0.01136 -0.00497 -0.01944 
H |g? -0.03222 -0.04345 =0.05324 -0.06169 -0.06887 
T RI= =: (116) -0.07487 -0.07975 -0.08358 -0.08640 -0.08828 
FRE -0.08928 -0.08944]; 
Xm = toeplitz(Xm11); % M-energy 
APPENDIX E Rr11 = 0.1*[1.77456 1.77298 1.76826 1.76042 ... 
1.74947 1.73548 1.71847 1.69854 1.67573 
MOM DATA 1.65016 1.62190 1.59106 1.55777 1.52213 
: : 1.48430 1.44439 1.40257 1.35897 1.31376 
MATLAB functions, scripts and data for planar rectan- 1.26710 1.21916 1.17009 1.12008 1.06929 
gles can be downloaded from http://www.eit.Ith.se/index.php? 1.01789 0.96607 0.91398 0.86181 0.80971 
0.75785 0.70639]; 


puid=175&projectpage=135&L=1. The zip-file contains the 


‘ Rr = toeplitz(Rril1)+eye (N) *5e-6; 
MATLAB functions F = eta0«(-lixkl)/4/pixones(1,N)*dx; % far field 


Calc_RX_matrices_rec(k1,k2,kN, 1x, ly, Nx, Ny, lib) B. Strip dipole with ly = 50£, = 0.48 and Ny = 100 
% calculates parts of the Xe,Xm,Rr matrices 
for a 1x,ly-rectangle and wavenumbers 

linspace (k1,k2,kN) 


ole 





% the data is stored in the folder lib K1=0.48*2*pi; k2=kl; kN=1; lx=1; ly=0.02; 
[X,p] = Nx=100; Ny=1; % \ellx=0.48\lambda 
r : = A re 2 1 

RX_rec_sym2full (X11,X12,X22,BxN, ByN, txN, tyN, pdef) lib = --\data ý 6 path to the data Liprary 

% calculates the Xe, Xm, Rr matrices from their parts be ee ae eee Cie k2,kN, 1x, ee Nya lib); 

F = farfieldmatrix(k,bas,meshp, evh, rvh) % load the data files 

% calculates the far-field matrix F for the Tosd eee A 
direction rvh and polatization evh /RowRX_rec_x=1_y=0p02_Nx=100_Ny=1_k=3p0159"')); 

F = sphmodematrix(k,bas,meshp, sphn) 


ol 


calculates the F-matrix for the spherical ... C. Strip dipole L= 0.1, Ny = 32 
modes sphn 


[a,GoQai,m,GoQrgap,adiff] = 
AntennaGoQ_NewtonlIt (a, Xe, Xm, F, gap0, adiff0,m0, abl) 





% Maximizes G/Q using Newton iterations for the ... eta0 = 299792458 x 4e-Txpi; % free space impedance 
dual problem kl = 0.1 x 2xpi; % wavenumber, 0.1 lambda 
Nx = 32; % number of elements 
and the scripts N = Nx-1; % number of unknowns 
dx = 1/Nx; % rectangle length 
dy = 0.02; % ractangle width 
maxGoQ_CVX_l_strip Xell =le3*[7.55508 -2.73908 -0.76452 -0.14218 
% calculates the Xe,Xm,Rr matrices for a strip -0.05383 -0.02634 -0.01488 -0.00924 -0.00614 
% and used CVX to maximize G/Q -0.00428 =0. 00311 =0:00233 -0.00179 -0.00141 
maxGoQ_CVX_2_strip -0.00112 -0.00091 -0.00075 -0.00063 -0.00053 
% calculates the Xe,Xm,Rr matrices for a strip ... -0.00045 -0.00038 -0.00033 -0.00029 -0.00025 
geometry for a range of wavenumbers and used ... -0.00022 -0.00020 -0.00018 -0.00016 -0.00014 
CVX or Newton iterations to maximize G/Q -0.00013 -0.00011]; 
maxGoQ_CVX_2_rec Xe = toeplitz(Xell); % E-energy 
% calculates the Xe,Xm,Rr matrices fora... Xm11 = [1.41378 0.76478 0.31782 0.20212 0.14925 
rectangle geometry for a range of 0.11844 0.09816 0.08376 0.07298 0.06460 . 
wavenumbers and used CVX or Newton 0.05789 0.05239 0.04779 0.04388 0.04052 
iterations to maximize G/Q 0.03760 0.03503 0.03275 0.03071 0.02887 
0.02721 0.02570 0.02431 0.02304 0.02186 
The length Zx is normalized to łą = 1m and wavenumber is 0.02078 0.01976 0.01882 0.01793 0.01710 
$ n : n 4 : s 0.01632]; 
in units of 1/¢, ie., the dimensionless quantity kf, is used. ym = tospiitz (mii); $ peenerdy 


Rr11 = le-3«[7.70515 7.70486 7.70397 7.70248 


7.70040 7.69773 7.69447 7.69061 7.68616 
A. Strip dipole 0 = 0.48A, N, = 32 7.68112 7.67549 7.66927 7.66246 7.65507 
7.64709 7.63852 7.62938 7.61965 7.60934 
7.59845 7.58699 7.57495 7.56234 7.54915 
eo : 7.53540 7.52109 7.50621 7.49076 7.47476 
eta0 = 299792458 » 4e-Txpi; % free space impedance 7.45821 7.44110]; 


oe 


kl = 0.48 * 2xpi; wavenumber, 0.48 


Rr = toeplitz(Rr11)+eye (N) *le-8; 


ae B nber GE cicrenta F = eta0+(-li*kl)/4/pixones (1,N)»xdx; % far field 
N = Nx-1; % number of unknowns 

dx = 1/Nx; % rectangle length D. Strip dipole with lx = 500, = 0.1, and Nx = 100 
dy = 0.02; % ractangle width 
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k1=0.1%2*pi; k2=k1; kN=1; 1lx=1; ly=0.02; Nx=100; 
Ny=1; % 1x=0.1\lambda 
lib = '..\data'; % path to the data library 


Calc_RX_matrices_rec(k1,k2,kN,1x,ly,Nx,Ny, lib); 
% load the data files 

load (strcat (lib, os 
'\RowRX_rec_x=1_y=0p02_Nx=100_Ny=1_k=0p62832')); 


E. Planar rectangle with l, = €,/2, Nx = 64 and Ny = 32 


The MoM data for a planar rectangle can be downloaded 
from http://www.eit.lth.se/index.php?puid=175 &projectpage= 
135&L=1. The zip-file contains the matrix elements associated 
with the first basis functions in the # and y directions. Run 
the matlab function 





RX_rec_sym2full 


to construct the Xe, Xm and R matrices. The matrices are 
computed using an equidistant mesh with Ny and N, elements 
in the @ and y-directions, respectively. The data is evaluated 
for planar plates with length ¿x and width @, using Nk 
equidistant samples of the wavenumbers k from the interval 
[ki, k2], ie, Linspace (k1,k2,Nk). The file 


load Xm_rec_x=1_y=0p02_Nx=50_Ny=1_k=0p02_3p2_50_6 





contains data for a plate with length ¢, = 1 and width 4, = 
0.02 meshed using Nx x Ny = 50 x 1 elements. It is the 6** 
wavenumber from the samples, i.e., 


kk = linspace(0.02,3.2,50); 
k = kk (6); 


The far-field matrix F in (27) for the direction 7 and 
polarization ê is computed with the MATLAB script 


evh = [1 0 0]; % unit vector in the x-direction 
rvh = [0 0 1]; % unit vector in the z-direction 
F = farfieldmatrix(k,bas,meshp,evh,rvh) ; 


The structures bas and meshp contain the basis functions 
and mesh used to discretize the structure. The corresponding 
F matrix for the spherical modes are computed using 


sphn = 6; % x-directed electric dipole 
F = sphmodematrix(kl,bas,meshp, sphn) ; 
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